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TECHNICAL  REPORT  SUMMARY 


One  purpose  of  this  investigation  has  been  to  regionalize 
the  upper  several  hundred  kilometers  of  the  mantle  under  the 
Arctic  region,  Siberia  and  the  Eurasian  continental  area  using 
seismic  surface  waves.  A  second  purpose  has  been  to  develop  and 
test  efficient  computer  techniques  for  the  computation  of  accurate 
theoretical  seismograms  for  hypothetical  sources  located  within 
Eurasia  and  which  make  use  of  the  results  of  the  regionalization 
part  of  this  study.  These  theoretical  seismograms  are  complete 
in  the  sense  that  they  contain  both  body  and  surface  waves;  they 
can  be  applied  directly  in  a  discrimination  program  by  comparing 
the  theoretical  seismograms,  computed  for  both  earthquakes  and 
underground  explosions,  with  the  actual  recorded  seismogram. 

Regionalization.  The  program  of  regionalization  involves 
the  measurement  of  phase  velocities  for  Rayleigh  waves  travers¬ 
ing  the  regions  under  investigation.  The  phase  velocity  curves 
were  obtained  with  the  single-station  phase  velocity  method, 
which  is  described  in  the  main  text  of  this  report.  Briefly, 
the  method  involves  selecting  records  of  earthquakes  from  the 
library  of  the  World  Wide  Standardized  Seismographic  Network 
(WWSSN)  for  events  which:  (1)  occurred  within,  or  on  the  peri¬ 

meter  of,  the  area  of  interest,  (2)  produced  good  long-period 
surface  wave  records  at  WWSSN  stations  for  which  the 
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epicenter-to-station  lines  lie  within  the  regions  being  in¬ 
vestigated,  and  (3)  generated  good,  large  recordings  at  a  suf¬ 
ficient  number  of  stations  to  ensure  an  accurate  fault  plane 
solution.  When  a  suitable  earthquake  was  found,  an  extensive 
data  processing  and  data  reduction  system  was  applied  to  trans¬ 
form  the  data  into  phase  velocity  curves  for  each  of  the 
selected  epicenter-to-station  lines. 

The  complete  set  of  phase  velocity  curves  (for  fundamental¬ 
mode  Rayleigh  waves  recorded  by  the  WWSSN  instruments),  which 
are  then  used  for  the  regionalization  of  the  Arctic  region, 

Siberia,  and  the  Eurasian  continental  area,  consists  of  about 
50  dispersion  curves.  Based  on  this  data  set,  we  have  regionalized 
the  area  of  interest  by  making  use  of  an  assumption  that  large 
parts  of  the  region  with  similar  geophysics  or  basement  geology 
and  age  will  have  similar  upper  mantle  structures.  A  similar 
assumption  was  made  with  success  in  the  regionalization  of  the 
Pacific  Ocean  basin  (Kausel,  Leeds  and  Knopoff,  1974;  Leeds, 
Knopoff  and  Kausel,  1974). 

The  inversion  has  proceded  by  the  use  of  both  linearized 
and  non-linear  inversions  procedures  with  the  mantle  properties 
in  each  of  the  geographic  regions  as  unknowns.  Our  regionalization 
has  made  use  of  maps  of  basement  geology  of  Eurasia  coupled  with 
a  postulate  of  similarity  of  mantle  cross-section  for  regions  of 
similar  basement  age,  made  previously  in  small-scale  regional 
studies  (Knopoff,  1972). 
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From  our  regionali zation  studies  on  the  Eurasian  continent, 
we  draw  the  following  conclusions:  First,  the  properties  of  the 
mantle  of  the  sediment-covered  Siberian  and  Russian  platforms 
are  highly  consistent  with  those  of  the  Baltic  and  Siberian 
shields  and  with  ancient  shields  observed  worldwide  (Knopoff, 

1972).  Second,  the  Tibetan  plateau  has  an  extremely  thick 
crust,  perhaps  as  great  as  75  km  from  surface  to  Moho.  Further¬ 
more,  the  upper  mantle  under  the  Tibetan  plateau  has  high  seismic 
velocities,  indicating  the  absence  of  partial  melting  to  relatively 
great  depths.  This  can  be  accounted  for  by  emplacement  of  the 
Indian  shield  under  Tibet  during  the  collision  of  the  Asian  and 
Indian  plates.  Third  the  Alpide  folded  belt  of  Iran  and  Turkey 
has  a  very  well-developed  low-velocity  channel  in  the  mantle, 
with  good  contrast  to  the  lid  above,  implying  the  presence  of  a 
zone  of  partial  melting  at  a  depth  of  about  90  to  100  km.  Fourth, 
the  Mongolia  -Sinkiang  geophysical  province  has  a  well-developed 
low-velocity  zone  in  the  upper  mantle  and  the  more-or-less  stable 
part  of  Eastern  China  also  has  a  low  velocity  channel  at  typical 
depths . 

The  principal  problem  in  the  inversion  has  been  the  con¬ 
struction  of  the  boundaries  to  the  geophysical  provinces,  for  which 
only  incomplete  information  is  found  in  the  literature.  A  full 
resolution  of  this  problem  is  still  in  the  future,  but  the 
regionalization  used  to  date,  in  which  Eurasia  is  divided  into  6 
rather  large  regions, is  statistically  consistent  with  the  data 
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and  has  about  the  correct  number  of  degrees  of  freedom  in  the 
model  parameterization. 

Theoretical  Seismograms.  The  second  phase  of  our  work  has 
been  the  development  and  testing  of  techniques  for  the  computation 
of  realistic  theoretical  seismograms  at  the  WWSSN  and  HGLP  in¬ 
stallations,  for  earthquake  and  explosive  sources  within  the  areas 
under  investigation.  For  SH  waves  excited  by  both  earthquakes 
and  explosions,  we  can  construct  complete  seismograms  down  to  a 
period  of  10  sec.  (and  occasionally  to  shorter  periods).  The 
difficulty  of  calculating  diffraction  effects  due  to  the  presence 
of  lateral  heterogeneity  remains  as  an  important  unsolved  problem; 
as  long  as  the  waves  cross  regional  boundaries  at  conditions  re¬ 
mote  from  grazing  incidence,  our  theoretical  seismograms  are 
probably  quite  accurate.  The  extension,  development  and 
optimization  of  these  algorithms  and  computer  programs  for 
Rayleigh,  or  P-SV  waves  on  a  spherical,  gravitating  earth  has 
also  been  one  of  the  main  efforts  under  this  contract.  The 
efficient  construction  of  accurate,  theoretical  SH  seismograms  for 
realistic  models  of  the  earth's  structure  using  multimode  methods 
is  based  on  effective  dispersion  computations,  attenuation 
calculations,  structural  transformations,  and  computations  of 
eigenfunction  characteristics,  effects  of  sphericity,  and 
point-source  response;  we  have  studied  all  these  ingredients  in 
detail . 


Our  work  with  theoretical  seismograms  has  concerned  the 
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multi-mode  surface  wave  phases  Lg ,  Sa  and  Sn;  references  to 
this  work  are  appended.  For  discrimination  studies,  it  is  also 
desirable  to  have  body-wave  phases  on  the  theoretical  seismograms; 
this  has  been  one  major  thrust  of  our  most  recent  work;  another 
has  been  the  improvement  in  efficiency,  control  of  accuracy  and 
increased  power  and  flexibility  of  the  computational  algorithms 
and  computer  programs.  We  have  reported  our  ability  to  generate 
theoretical  seismograms  with  body  and  surface  waves, on  the  same 
record  using  up  to  twenty-one  modes.  Multimode  theoretical 
seismograms  containing  body  wave  phases  have  been  generated  in 
the  past  (Sato,  Usami  and  Landisman,  1968) ;  however,  these  time 
series  were  limited  to  ultralong  periods.  The  important  point 
of  our  recent  results  is  that,  owing  to  the  efficiency  of  our 
new  algorithms,  we  have  been  successful  in  extending  the  period 
content  of  the  theoretical  seismograms  through  the  range  covered 
by  the  WWSSNLP  instruments.  Thus,  we  can  generate  the  theoretical 
time  series  which,  for  the  first  time  we  believe,  permits  us  to 
compare  theoretical  seismograms  directly  with  the  entire  records 
obtained  at  the  WWSSN  and  the  HGLP  installations  down  to  a  period 
of  ten  seconds.  This  requires  that  90-100  modes  be  used  in  the  multi- 
mode  synthesis.  The  necessary  frequency-domain  information  is 
obtained  in  a  single,  relatively  inexpensive  computer  run;  time 
series,  for  any  source  specification,  are  then  obtained  with  a 
single  run  of  a  second  program  (Liao,  Schwab  and  Mantovani  1977; 
preprint  appended) . 
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In  the  case  of  Rayleigh  waves,  sphericity  and  gravity  intro¬ 
duce  complications  not  encountered  with  Love  waves.  For  Love 
waves,  optimized  flat-structure  programs  can  be  used  since  the 
sphericity  of  the  real  earth  can  be  handled  by  exact  transformation, 
and  since  gravity  does  not  affect  Love  waves.  At  present,  accurate 
Rayleigh-wave  computations  for  a  real  earth  require  that  sphericity 
and  gravity  be  handled  directly  in  the  computational  algorithms 
and  programming;  hence,  our  first  task  here  was  to  improve  the 
accuracy  characteristics  of  direct  Rayleigh-,  or  spheroiaal-wave 
computations.  The  difficulty,  expense,  and  time  required  for  this 
analysis  showed  us  why  this  information  was  not  in  the  literature, 
even  though  the  basic  algorithm  is  available  (Alterman,  Jarosch 
and  Pekeris,  1959).  The  results  of  our  work  (Schwab  et  al ,  1977) 
are  appended.  They  include:  (1)  an  improved  and  simplified  com¬ 

putational  algorithm  for  the  computation  of  the  phase  velocity  of  Rayleigh 
waves.  (2)  What  we  believe  to  be  the  first  direct  method  for 
computing  the  group  velocity,  i.e.,  without  appeal  to  variational 
methods;  (3)  numerical  analyses  and  examples  of  numerical  dif¬ 
ficulties  encountered  in  this  type  of  computation;  (4)  detailed 
analysis  of  the  efficiency  of  our  optimized  algorithm  and  pro¬ 
gramming.  Even  in  this  optimized  form,  these  direct  Rayleigh- 
wave  computations  for  a  spherical,  gravitating  earth,  are  about 
six  times  slower  than  the  comparable  Love-wave  computations; 
whereas  we  know  that  computations  for  Rayleigh  waves  on  a  flat, 
non-gravitating  structure  are  only  twice  as  slow.  Our  conclusion, 
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from  the  Rayleigh  wave  work  done  here,  is  that  new  algorithms 
are  required  before  the  efficiency  of  the  Rayleigh-wave  com¬ 
putations  (including  sphericity  and  gravity)  will  approach 
a  level  justifying  computation  of  the  associated  multimode, 
theoretical  seismograms  down  to  the  desired  period  of  ten 
seconds.  We  have  developed  a  sufficiently  improved  algorithm 
(Schwab,  1977) ,  but  the  necessary  numerical  verification  and 
testing  are  only  in  their  initial  stages. 

References  to  our  own  work  on  theoretical  seismograms  for 
Sa,  Lg  and  Sn  include:  Knopoff,  Schwab  and  Kausel  (1973); 
Knopoff  et  al  (1974);  Nakanishi,  Schwab  and  Kausel  (1976) ; 
Kausel,  Schwab  and  Mantovani  (1977);  Mantovani  et  al  (1977). 
References  to  our  work  on  theoretical  seismograms  which  in¬ 
clude  both  body  and  surface  waves  are:  Nakanishi,  Schwab 
and  Kausel  (1976) ;  Kausel,  Schwab  and  Mantovani  (1977); 
Mantovani  et  al  (1976);  Mantovani  (1977a, b);  Liao,  Schwab 
and  Mantovani  (1977). 
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TECHNICAL  REPORT 


I.  Purposes 

'One  purpose  of  this  project  was  to  determine  the  regional  variations 
of  the  crust  and  upper  mantle  in  the  regions  under  investigation 
using  the  properties  of  surface  wave  dispersion  in  the  Arctic 
region  and  the  Eurasian  continental  area.  The  regions  of  high 
seismicity  within  and  around  the  region, plus  the  dense  set  of 
WWSSN  stations  located  around  the  perimeter  of  the  region,  gave 
assurance  that  we  could  acquire  sufficient  single-station  phase 
velocity  data  for  the  application  of  regionalization  procedures. 

A  second  purpose  was  to  develop  computational  techniques 
for  calculating  complete,  accurate  theoretical  seismograms  for 
both  earthquake  and  explosion  sources  within  the  area  under  in¬ 
vestigation,  and  which  could  be  compared  with  those  recorded  at 
the  WWSSN  and  HGLP  stations  at  the  edge  of  this  region. 

II.  Review  of  scientific  background 

Regionalization.  The  single-station  surface  wave  method 
allows  all  stations  to  be  located  at  the  edge  of  the  region 
under  investigation;  this  technique  is  ideal  for  the  study  of 
Eurasia.  Knopoff  and  Schwab  (1968)  extended  the  description  of 
the  single  station  method  (Brune,  Nafe  and  Oliver,  1960)  to  take  into 
account  the  frequency  dependence  of  the  apparent  initial  phase 
of  the  source.  The  success  of  the  single-station  method 
depends  on  the  calculation  of  the  initial  phase.  For  dislocation 
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sources  in  layered  media,  the  initial  phase  can  be  determined 
by  a  method  due  to  Ben  Menahem  and  Harkrider  (1964)  and 
Harkrider  (1964,  1970)  who  obtained  the  surface  wave  response  to 
double  couple  sources.  The  far-field  response  to  displace¬ 
ment-dislocation  faulting  is  equivalent  to  that  from  a  point- 
source,  double-couple  in  an  unfaulted  medium  (Burridge  and 
Knopoff,  1964).  Thus  if  the  source  focal  mechanism  is  known 
from  fault-plane,  surface  wave  amplitudes,  or  other  methods, 
the  initial  phase  can  be  determined  for  small  or  moderate 
sized  earthquakes.  By  means  of  transformation  techniques 
(Biswas  and  Knopoff,  1970)  or  empirical  correction  methods 
(Bolt  and  Dorman,  1961)  ,  it  is  possible  to  provide  corrections 
for  sphericity. 

The  first  inversions  of  surface  wave  dispersion  data  to 
give  upper  mantle  structure  were  carried  out  by  Dorman  and 
Ewing  (1962),  Brune  and  Dorman  (1964)  and  Knopoff,  Mueller  and 
Pilant  (1966)  .  These  papers  were  concerned  with  obtaining  a 
single  structure  which  fit  the  experimental  data.  Knopoff 
(1961,  1962)  pointed  out  that  the  inversion  of  noise-free 
data  is  not  unique.  Subsequent  efforts  were  mainly  concerned 
with  finding  the  set  of  structural  models  which  fit  the  data  to 
within  the  experimental  accuracy  (Keilis-Borok  and  Yanovskaya, 
1967;  Press,  1968,  1969).  The  inversion  of  noisy  data  enlarges 
the  span  of  non-unique,  acceptable  models  over  that  for  noise- 
free  data. 

Two  methods  of  inversion  are  available.  In  the  first  of  the 
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techniques,  perturbation  methods  are  applied  in  a  linearized 
version  of  the  problem.  The  problems  of  inconsistent  or 
correlated  data  are  accounted  for  by  smoothing  in  the  data- 
space,  usually  by  a  least-squares  procedure.  The  problems  of  over¬ 
parameterization  of  the  model  space  are  taken  care  of  by  smoothing 
in  the  model  space,  by  a  procedure  usually  described  as  "deltaness". 
Further  reductions  of  the  number  of  model  parameters  can  stabilize 
the  inversions.  Examples  of  linearized  inversions  are 

given  by  Backus  and  Gilbert  (1968,  1970),  Knopoff  &  Jackson  (1972) 
Jackson  (1972),  and  Wiggins  (1972). 

The  full  non-linear  inverse  problem  in  a  multidimensional 
parameter  space  of  high  order  has  been  attacked  by  Monte  Carlo 
methods  (Press,  1968,  1969).  The  inefficiency  of  Monte  Carlo 
methods  can  be  minimized  at  the  sacrifice  of  reducing  the 
dimensionality  of  the  parameterization  to  a  value  less  than  about 
seven  by  exploring  a  neighborhood  of  acceptable  solutions  for  other 
acceptable  solutions.  Such  a  program  has  been  called  Hedgehog 
and  has  been  much  used  in  this  laboratory  (e.g.  Biswas  and 
Knopoff,  1974  ;  Knopoff  and  Schlue,  1972  ;  etc.);  copies  of 
this  program  exist  in  Moscow,  Bologna ,  Edmonton,  Bari,  Cambridge  (IK),  Paris, 
etc.  Both  linearized  and  Hedgehog  methods  have  been  used  in  our 
inversions . 


11 


The  single-station,  surface  wave  regionalization  of  an 
area  is  based  on  the  assumption  that  the  surface  waves  travel 
paths,  and  the  assumption  that  the  phase  travel  time  from  epi¬ 
center  to  station  is  the  sum  of  the  travel  times  through  the 
homogeneous  subdivisions  of  the  laterally  heterogeneous  region 
(Knopoff,  1969).  This  permits  the  travel-time  at  each  frequency 
to  be  expressed  as  a  system  of  inhomogeneous  equations; 
for  each  path,  the  total  travel  time  is  the  inhomogeneous  term, 
the  distances  through  the  subdivisions  are  the  coefficients,  and 
the  slownesses  in  the  subdivisions  are  the  unknowns.  Ke  have 
pointed  out  (Leeds  et  al .  1974)  that  the  solution  of  this  system 
of  equations  can  yield  the  experimental  phase  slownesses 
associated  with  each  of  the  subdivisions  only  if  we  make  the 
assumption  that  the  errors  in  each  region  are  uncorrelated. 

Since  these  errors  are  not  uncorrelated,  we  must  consider  the 
model  parameters  as  the  primary  unknowns  in  the  inversion  and 
derive  the  phase  slownesses  in  the  pure  regions  therefrom. 

The  first  application  of  the  use  of  single-station  phase 
delay  measurements  in  the  pursuit  of  regionalization  of  a  large 
inhomogeneous  area  was  performed  by  Kausel,  Leeds  and  Knopoff 
(1974),  Leeds,  Kausel  and  Knopoff  (1974)  and  Leeds  (1975).  In 
the  above  work,  phase  delays  along  long  paths  across  the  Pacific 
area  were  found  to  vary  systematically  with  distance  from  the 
East  Pacific  Rise.  This  persuaded  us  that  the  appropriate 
geographic  regionalization  was  plausibly  based  on  a  basement 
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i.e.  magnetic,  age  of  sea-floor  spreading.  Our  geographic 
provinces  then  took  the  form  of  broad  magnetic  age  stripes. 

An  inversion  using  this  regionalization  showed  that  although  the 
number  of  criss-crossing  paths  was  plentiful,  the  number  of 
degrees  of  freedom  in  the  data  was  remarkably  small.  furthermore, 
we  found  that  the  data  set  did  not  permit  one  to  obtain  detailed 
information  concerning  the  bottom  of  the  low-velocity  channel.  We  did 
find  that  the  lithosphere  increased  in  thickness  monotonically  with  age: 
the  ridge  crest,  the  lid  of  the  channel  has  almost  zero  thickness, 
while  in  the  oldest  ocean  this  lid  is  about  100  km  thick.  The  lid  thick 
nesses  are  consistent  with  a  geochemical  model  in  which  the  lid- 
channel  interface  is  at  the  solidus  for  wet  peridotite. 

Because  they  have  not  developed  good  long-period  instruments, 

Soviet  seismologists  have  not  been  able  to  focus  attention  on 
upper  mantle  studies.  Also,  their  extensive  program  in  deep  seismic 
sounding  has  focused  interest  upon  crustal  studies.  The  Soviet 
surface  wave  work  has  been  limited  to  short-period  investigations  — 
usually  less  than  40  seconds  --  and  has  been  concerned  mainly  with 
determining  crustal  properties.  References  to  Soviet  surface-wave 
work  include  Arkhangelskaya  (1960),  Savarensky  and  Ragimov  (1958, 

1959),  Savarensky,  Solov'eva  and  Shechkov  (1959)  and  Savarensky 
and  Sikharulidze  (1959),  Popov  (1960),  Shechkov  (1961,  1964,  1970), 
Savarensky  and  Shechkov  (1961),  Shechkov  and  Solov'eva  (1961), 

Savarensky  and  Peshkov,  1968;  Sikharulidze  and  Makharadze,  1968; 
Savarensky  et  al ■  1969. 
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Since  appropriate  Soviet-b^sed  seismic  data  are  not  avail¬ 
able  to  us,  our  own  work  has  thus  far  focused  on  the  estimation 
of  upper  mantle  properties  in  Eurasia  through  the  use  of  the 
only  reasonable  tools  available  to  us,  namely  phase  delays  of 
surface  waves  on  long  paths,  criss-crossing  the  region  of 
interest.  This, plus  a  plausible  model  of  regionalization  -- 
based  in  part  on  the  observations  in  Knopoff  (1972)  for 
continents  — .permits  an  attack  on  the  problem  of  regionalization 
of  the  upper  mantle  of  Eurasia. 

Theoretical  seismograms.  Relative  to  the  discrimination 
problem,  probably  the  most  important  feature  in  the  calculation 
of  theoretical  seismograms  which  requires  improvement  over 
previously  existing  systems  for  such  computations  is  the  capa¬ 
bility  of  extending  both  body-  and  surface-wave  portions  of  the 
theoretical  computer  seismograms  to  short  periods.  In  this 
context,  by  "short-period"  we  refer  to  the  period  range  covered 
by  the  long-period  instruments  of  the  WWSSN  installations.  In 
our  formulation,  the  successful  accomplishment  of  this  task  is 
dependent  upon  improved  techniques  for  obtaining  multimode 
dispersion-attenuation  information  for  reasonably  realistic 
models  of  the  earth,  i.e.  spherical,  radially  heterogeneous, 
anelastic  models. 

Optimization  of  this  type  has  been  one  of  the  main  interests 
in  our  laboratory  for  several  years.  The  results  of  our  early 
work,  based  on  the  Thomson  ( 1950) -Haskell  (1953)  technique  and 
on  Knopoff' s  (1964)  method  for  treating  flat-layered  structures. 


are  reported  by  Schwab  (1970)  and  Schwab  and  Knopoff 
(1970;  1971;  1972;  1973).  The  results  of  our  work  on 
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spherical-to-f lat-structure  transformation  techniques,  which 
permit  the  use  of  flat-structure  programs  in  dispersion- 
attenuation  computations  with  Love  waves  for  spherical  models 
of  the  earth,  are  given  by  Biswas  and  Knopoff  (1970),  Schwab 
and  Knopoff  (1971;  1972;  1973),  and  Kausel  and  Schwab  (1973). 

In  this  last  reference,  we  have  also  given  an  outline  of  the 
approach  we  have  adopted  to  handle  the  synthesis  of  multimode 
seismograms  once  the  dispersion,  attenuation,  source,  and 
excitation  functions  have  been  specified.  The  entire 
theoretical  seismogram  for  a  dislocation  source  in  a  spherical 
earth  can  be  expressed  as  a  simple  sum  of  normal  mode  contri¬ 
butions  (Saito,1967;  Takeuchi  and  Saito,  1972).  We  first  applied 
our  scheme  for  generating  theoretical  seismograms  to  the  inter¬ 
pretation  of  the  seismic  phase  Lg  (Knopoff,  Schwab  and  Kausel, 

1973;  Knopoff  et  al .  1974)  This  phase  is  a  multimode  interference 
pheonomenon  which  belongs  to  the  surface  wave  portion  of  the  seismo¬ 
gram.  In  this  report  we  indicate  that  we  have  developed  a  system 
for  synthesis  where  both  body  and  surface  waves  appear  on  the 
same  seismogram  for  realistic  models  for  the  earth,  and  where  the 
period  range  spans  that  covered  by  the  long-period 
instruments  of  the  WWSSN.  Earlier  work  of  this  type,  which  was 
performed  with  simplified  models  of  the  earth,  is  summarized  by 
Alterman  and  Loewenthal  (1972).  Sat6,  Usami  and  Landisman  (1968) 
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describe  the  computation  of  com;  lete  theoretical  seismograms 
for  realistic  models  of  the  earth.  However,  their  results  are 
limited  to  ultralong  periods. 

III.  Objectives,  methods  and  results. 

Data  collection.  Accurate  surface  wave  phase  delays 
across  the  area  with  single-station  method  require  that  we 
know  the  focal  mechanism  and  the  depth  of  focus  in  order  to 
obtain  the  corrections  due  to  the  apparent  initial  phase.  Vie 
have  confined  our  attention  to  the  measurement  of  the  phase 
velocity  dispersion  of  the  fundamental  mode .  I'rez  and  Schwab 
(1976)  have  computed  the  importance  of  the  structural  parameters 
on  the  determination  of  the  initial  phase. 

Long  period  records  from  the  47  WWSSN  stations  which 
border  the  region  of  interest  have  been  used  in  the  study. 

The  locations  of  these  stations  are  shown  in  Figure  1.  We  have 
also  obtained  several  seismograms  (through  World  Data  Center  B) 
of  Soviet  records  made  in  Central  Asia  on  experimental  long- 
period  instruments.  However,  the  precision  of  these  latter 
recordings  plus  an  uncertain  calibration  impulse  response  has 
not  permitted  us  to  use  them. 

An  example  of  the  intermediate  magnitude  seismicity  of 
the  area  is  given  in  Figure  2.  Epicenters  are  plotted  for 
earthquakes  which  occurred  during  the  interval  from 
February,  1963  to  February,  1967  and  having  magnitudes  between 
5.9  and  6.6.  In  addition  to  the  epicenters  shown  in  Figure  2, 
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FIG.  1.  Locations  of  WWSSN  and  HGLP  stations. 

FIG.  2.  Seismicity  of  the  Eurasian  region.  The 

regions  of  high  seismicity  along  the  eastern 
border  of  the  Kamchatka  peninsula  and  along 
the  Aleutian  arc  also  provide  useful  events 
for  the  study. 
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there  are  regions  of  high  seismicity  along  the  eastern  border 
of  the  Kamchatka  peninsula  and  along  the  Aleutian  arc.  The 
choice  of  this  range  of  magnitudes  is  governed  by  two  considerations. 
First,  experience  has  shown  that  good  long-period  surface  wave 
information  requires  events  with  magnitudes  above  a  certain  value; 
of  course,  a  shock  which  is  so  large  as  to  send  the  instrument 
off-scale  is  useless  for  our  purposes.  Second,  the  application 
of  the  single-stations  method  requires  knowledge  of  the  focal 
mechanism.  We  must  therefore  use  events  large  enough  to  allow 
us  to  obtain  an  accurate  fault  plane  solution  for  each  event  we 
select  for  processing. 

Since  the  set  of  stations  around  the  area  to  be  studied 
is  dense  as  are  the  regions  of  high,  intermediate-magnitude 
seismicity  located  within  and  around  the  area,  there  has  been 
no  problem  in  obtaining  sufficient  data  for  the  project.  It  is 
interesting  to  note  that  the  area  is  almost  completely  encircled 
by  either  stations  or  epicenters  or  both.  The  limits  of  the 
area,  which  we  have  covered  with  a  sufficiently  dense  set  of 
paths  from  earthquakes  to  epicenters,  are  shown  in  Figure  3 
by  the  solid  lines.  The  shaded  regions  are  those  of  high,  large- 
magnitude  seismicity. 

Recordings  from  fifteen  events  were  processed.  The  list 
of  events  is  given  in  Table  I,  with  their  USCGS-NOAA  specifi¬ 
cations.  The  focal  parameters  and  other  data  are  listed  in 
Table  II.  We  have  constructed  fault  plane  solutions  for  11 
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FIG.  3.  Limits  (solid  lines)  of  the  region  covered 

with  a  dense  set  of  epicenter-to-station  paths. 
Solid  regions  are  those  of  high  seismicity. 


Table  I 


Events  used  in  phase  velocity  determinations 
(USCGS  specifications) 


1.  Red  Sea 

2.  Hsingtai 

3.  Kamchatka 

4 .  Lena  River 
'  Tashkent 

6 .  Lop  Nor 

7.  Yunnan 

8.  Hindu  Kush 

9.  Hindu  Kush  : 

10.  East. 


Time 

Lat. 

Long. 

h  (km) 

M 

,  1969 

07:15:54.4 

27 . 7°N 

34 . 0°E 

33 

6.0  to 
6.8 

,  1966 

21:29:17.4 

37 . 3°N 

114. 9°E 

33 

6.0 

,  1964 

14:30:29.1 

51 . 8°N 

156. 8°E 

136 

5.7 

,  1964 

13:47:20.6 

78 . 2°N 

126. 6°E 

50 

6.1 

■  1966 

23:22:49.3 

41 . 3°N 

69 . 2°E 

8 

5.  3 

1970 

07:29:58.6 

40 . 9°N 

8  9 . 4  °E 

0 

1966 

10:44:41.3 

26 . 1°N 

103. 2°E 

33 

5.7 

1966 

07:46:16.1 

36 . 4°N 

71 . 1°E 

221 

6.2 

1974 

12:11:43.7 

35 . 1°N 

72 . 9°E 

22 

6.0 

1965 

01:40:33.2 

53 . 2°N 

161. 9°W 

33 

6.4  to 
6.7 

1965 

16:50:23.6 

53 . 3°N 

161. 8°W 

33 

6.1  to 
6.6 

1965 

09:32:09.3 

52 . 3°N 

174. 3°E 

41 

5.9  to 
6.5 

1965 

04:02:53 

52 . 1°N 

175. 7°E 

35 

5.9  to 
6.2 

18 


Table  II 

Focal  Parameters  for  Events  Studied 


Event  No. 

$ 

5 

h  (km) 

Fault  Plane 
Solutions 

Pa  this 

1. 

133° 

66° 

242° 

15 

Fig.  4 

Fig.  14 

2. 

122° 

82° 

0° 

14.5 

5 

15 

3. 

210° 

90° 

90° 

136. 

6 

16 

4  . 

165° 

5  8°W 

260oa,b 

11 . 5a ' b 

7 

17 

5. 

305°d 

7  0od 

90od 

8. 

22 

6  . 

explosion 

0. 

22 

7. 

190° 

0 

o 

279° 

33. 

8 

18 

8  . 

77° 

cn 

o 

0 

76° 

221. 

9 

18 

9  . 

122° 

0 

o 

■^r 

61° 

22. 

10 

18 

10. -11.  (twin  earthquakes)  11a, b  19 

12.-13.  (twin  earthquakes)  12  20 

a.  Depth  obtained  from  Rayleigh  wave  spectra. 

b.  Rayleigh  wave  spectra  for  event  2  shown  in  Fig.  13 

c.  Sykes  (1967)  gives  fault  planes  for  this  event  as 

v  =  4°,  6  =  58°W 

4  =  338°, 


6  =  54°E 
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of  these  cases  (including  two  pairs  of  twin  earthquakes) .  The 
list  of  figures  of  these  fault  plane  solutions  is  given  in 
Table  II.  In  one  case  (No.  5)  we  have  used  the  fault  plane 
solution  of  Zakharova  et  al  (1971).  In  one  case  (No.  6)  we 
have  assumed  the  source  that  of  an  explosion.  In  two  cases 
(Nos.  1,  2)  we  have  used  Rayleigh  wave  spectra  to  help  deter¬ 
mine  the  focal  parameters;  in  these  cases,  we  were  able  to 
obtain  the  depth  of  focus  more  accurately  than  by  conventional 
methods.  In  the  case  of  No.  2  the  spectra  (Fig.  13)  also  helped 
refine  the  fault  plane  parameters  and  permit  us  to  modify  Sykes' 
(1967)  values.  The  methods  used  to  determine  fault-plane 
parameters  from  surface  wave  spectra  are  given  by  Ben  Menahem 
and  Toksoz  (1963). 
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FIG.  4.  Fault  plane  solution  for  event  4  (Red  Sea) 

occurring  at  07:15:54.4  GMT,  March  31,  1969. 

FIG.  5.  Fault  plane  solution  for  event  1  (Hsingtai) 
occuring  at  21:29:17.4  GMT,  March  7,  1966. 

FIG.  6.  Fault  plane  solution  for  event  3  (Kamchatka) 

occuring  at  14:30:29.1  GMT,  December  26,  1964. 

FIG.  7.  Fault-plane  solution  for  event  2  (Lena  River) 
occurring  at  13:47:20.6  GMT,  August  25,  1964. 

The  solution  given  by  Sykes  (1967)  is  indicated 
by  dotted  lines. 

FIG.  8.  Fault  plane  solution  for  event  7  (Yunnan) 

occurring  at  10:44:41.3  GMT,  February  13,  1966. 

FIG.  9.  Fault  plane  solution  for  event  8  (Hindu  Kush) 
occurring  at  07:46:16.1  GMT,  June  6,  1966. 

FIG.  10.  Fault  plane  solution  for  event  9  (Hindu  Kush) 
occurring  at  12:11:43.7  GMT,  December  28,  1974. 

FIG.  11a.  Fault  plane  solutions  for  events  10  and  11 
''■''■k'  (Eastern  Aleutians)  occurring  at  01:  40:33.2 

(Fig.  11a)  and  16:50:29  GMT  (Fig.  11b),  February 
1965. 
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FIG.  12.  Fault  plane  solution  for  events  12  and  13 

(Western  Aleutians)  occurring  at  09:32:09.3  GMT, 
February  5,  1965  and  04:02:53  GMT,  February  6, 
1965. 

FIG.  13.  Rayleigh-wave  amplitude  distributions  for 

event  2  (Lena  River)  occurring  at  13:47:20.6 
GMT,  August  25,  1964.  The  central  set  of 
radiation  patterns  are  the  results  of  theor¬ 
etical  computations  based  on  the  fault  plane 
solution  given  in  the  text.  The  other  four 
radiation  patterns  depict  the  experimental 
»  results. 
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Two  points  should  be  noted  concerning  the  accuracy  of  the 
digitizations  of  the  recorded  events.  First,  most  of  the  event 
records  we  have  used  are  about  as  large  as  they  could  be  without 
going  off  scale.  This  has  necessitated  a  change  in  our  data 
processing  techniques  which  should  be  noted  for  the  information 
of  others  involved  in  this  type  of  work. 

In  the  past,  our  standard  procedure,  when  working  with 
smailer-amplitude  recordings,  has  been  to  digitize  from  copies 
of  35  mm  microfilm  records  of  the  WWSSN  seismograms  made  with  a 
standard  microfilm  reader-printer  (Itek  18.24  Reader-Printer). 
Tests  which  compare  the  phase  velocity  results  obtained  from 
full-size  record  copies  provided  by  NOAA  with  the  results 
obtained  from  our  microfilm  copies  show  that  distortion  in  the 
copying  process  is  of  concern  when  working  with  large-amplitude 
recordings  such  as  those  employed  in  the  present  study.  We  en¬ 
courage  the  use  of  full-size  record  copies  of  large  events 
obtained  directly  from  NOAA. 

The  second  point  concerns  the  fact  that  the  direction  of 
swing  of  the  galvanometer  may  not  be  parallel  to  the  axis  of 
the  recording  drum.  Although  Jaimes  and  Linde  (1971)  term  this 
phenomenon  "a  source  of  major  error  in  digital  analysis  of  WWSSN 
seismograms",  our  tests  show  the  effect  to  be  negligible  on  phase 
travel  times  computed  using  the  single-station  method  for  epicenter- 
station  separations  of  a  few  thousand  kilometers.  In  the  case  of 
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the  poorest  galvanometer  alignment  we  encountered,  about  triple 
the  normal  ramp  slope  of  0.3°,  we  found  only  negligible  differences 
between  the  phase  velocity  curve  obtained  with  the  correct 
digitization  base  line,  and  the  curve  obtained  with  the  normal 
ramp  as  a  base  line. 

The  paths  over  which  the  phase  delays  have  oeen  measured  from 
each  of  the  seismic  events  are  shown  in  Figs.  14-20  and  are  sum¬ 
marized  in  Table  II.  Sample  phase  velocity  results  are  given  in 
Figure  21,  which  illustrates  the  variation  in  dispersion  for  dif¬ 
ferent  propagation  paths.  V.’e  have  obtained  phase  velocity  data  for 
most  paths  over  a  period  range  extending  from  about  30  or  38  sec. 
in  most  cases,  to  as  long  as  357  sec  in  a  few  rare  cases.  The 
instrumental  response  at  these  longest  periods  is  unreliable;  if 
the  measured  values  of  phase  velocity  and  its  first  derivative 
were  in  significant  disagreement  with  values  from  the  free  mode 
spectrum,  then  these  long  period  values  were  rejected.  For  our 
present  inversions,  we  have  only  used  periods  up  to  250  sec.  The 
specific  period  ranges,  which  we  have  used  in  the  inversions  are 
shown  in  Table  3. 

The  phase  velocities  from  five  earthquakes  and  one  nuclear 
explosion  for  the  32  paths  crossing  Eurasia  (Figure  22)  sort 
themselves  into  two  groups  (Figures  23  and  24).  The  paths  with 
higher  phase  velocities  are  generally  those  that  cross  the  stable 
platforms  and  shields  (such  as  paths  from  the  Hsingtai  earthquake 
to  Scandinavian  (KEV)  and  German  (STU)  stations;  typical  of  the 
lower-velocity  group  is  the  phase  velocity  on  the  paths  from  the 
Red  Sea  earthquake  to  the  southern  Asiastic  stations  (SHL,  MAN,  etc.)). 
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FIG. 


FIG. 


FIG. 


FIG. 


FIG. 


14.  Epicenter-to-station  paths  processed  for 
event  4  (Red  Sea)  occurring  at  07:15:54.4  GMT, 
March  31,  1969. 

15.  Epicenter-to-station  paths  processed  for 
event  1  (Hsingtai)  occurring  at  21:29:74.4  GMT, 
March  7,  1966. 

16.  Epicenter-to-station  paths  processed  for 

event  3  (Kamchatka)  occurring  at  14:30:29.1  GMT, 
December  26,  1964. 

17.  Epicenter-  to-station  paths  processed  for 

event  2  (Lena  River)  occurring  at  13:47:20.6  GMT, 
August  25,  1964. 

18.  Epicenter-to-station  paths  processed  for  event 
7  (Yunnan)  occurring  at  10:44:41.3  GMT, 

February  13,  1966,  and  events  8  and  9 

(Hindu  Kush)  occurring  at  07:46:16.1  GMT, 

June  6,  1966  and  12:11:43.7  GMT,  December  28, 


1974. 
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FIG. 


FIG. 


FIG. 


19.  Epicenter-to-stacion  paths  processed  for 
events  10  and  11  (Eastern  Aleutians)  occurring 
at  01:40:33.2  (dashed  lines)  and  16:50.29  GMT 
(solid  lines),  February  6,  1965. 

20.  Epicenter-to-station  paths  processed  for 
events  12  and  13  (Western  Aleutians)  occurring 
at  09:32:09.3  GMT,  February  5,  1965  (dashed 
line)  and  at  04:02:53  GMT,  February  6,  1965 
(solid  line) . 

21.  Sample  phase  velocity  results  for  the  paths 
from  the  Hsingtai  earthquake  (1966)  to  STU 
and  from  the  Lena  River  earthquake  (1964)  to 
HOW.  The  paths  are  shown  in  Figures  15  and 
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TABLE  3 

PERIOD  RANGE  (sec) 


PATH 

VUNAN-ANP 

-BAG 

-MSH 

-SEO 

-SHK 


30  38  50  69  100  119  139  167  192  208  227 

X _ X 

X _ X 

X _ X 

X  X 

X  X 


HINDU  KUSH  I-ANP 
-HKC 


X 

X 


X 


X 


HINDU  KUSH  I I-ANP 


X 


X 


1 


250 
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FIG.  22.  Propagation  paths  across  the  Eurasian  continent 
from  five  earthquakes  and  one  nuclear  explosion 

FIG.  23.  All  Eurasian  phase  velocities  can  be  sorted  into 
two  groups  (shaded  areas) ,  except  for  phase 
velocities  POO-2  and  SHL-1,  which  fall  between 
these  two  groups.  Phase  velocities  for 
"standard"  shield  (FLO-GOL) ,  younger  stable 
regions  (SHA-LUB)  and  rift  zones  (TUC-BOZ) 
are  shown  for  comparison  (Biswas  and  Knopoff, 
1974).  The  global  average  phase  velocities 
obtained  from  free-mode  observations  are  also 
shown  (F.M.)  (Gilbert  and  Dziewonski,  1975). 

FIG.  24.  Propagation  paths  corresponding  to  the  two 

phase  velocity  groups  in  Fig.  21.  The  solid 
line  indicates  a  path  with  higher  phase  velocity; 
the  dashed  line  indicates  a  path  with  lower 
phase  velocity. 
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Phase  velocities  for  typical  shield  regions  and  young  stable  con¬ 
tinental  regions  (Biswas  and  Knopoff,  1974)  are  shown  for  com¬ 
parison.  The  incompatibility  of  most  of  the  phase  velocity 
observations  for  Eurasian  paths  with  the  "standard"  curves  for 
homogeneous  regions  testifies  to  the  inhomogeneous  nature  of  the 
Eurasian  region. 

Regionalization 

The  first  step  in  the  structural  analysis  requires  the 
division  of  the  area  into  subregions;  each  of  these  is  assumed 
to  be  laterally  homogeneous.  In  our  first  attempt  we  have  sub¬ 
divided  the  Eurasian  area  into  six  broad  regions  (Fig.  25) .  These  are 


1. 

Ancient  PreCambrian 

Shields 

2. 

Stable  Platforms 

3. 

The 

Himalayan-Alpide 

Mountain  Belt 

4. 

The 

Tibetan  Plateau 

5. 

The 

Sinkiang-Mongolian  Seismic  Zone 

6  . 

The 

Chinese  "Stable" 

Region 

This  choice  of  regionalization  is  largely  governed  by  an  attempt 
to  define  a  small  number  of  discrete  geographically  homogeneous 
provinces.  In  this  choice  we  have  been  guided  by  tectonic  maps 
of  Eurasia  (Khain  and  Muratov,  1969)  as  well  as  by  seismicity  and 
sparse  heat  flow  information (Lubimova  and  Polyak,  1969).  Our  first 
impulse  was  to  regionalize  the  stable  parts  of  Eurasia  according  to 
the  bimodal  division  of  these  regions  into  stable  ancient  shields 
and  stable  younger  continents  (Knopoff,  1972).  However, 
information  about  basement  ages 
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FIG. 


25.  Regionalization  of  the  Eurasian  continent 

based  on  the  tectonic  map  of  Khain  and  Muratov 
(1969)  . 
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of  the  Eurasian  sediment-covered  platforms  is  difficult  to  obtain, 
so  we  have  taken  the  simple  expediment  of  dividing  the  stable 
parts  of  the  continent  into  the  geologically  identif icable  shields 
and  a  second  region  which  includes  the  rest  of  the  seisniically  stable  land 
mass;  these  form  our  regions  1  and  2.  In  selecting  region  2  as  a 
homogeneous  region  we  impose  no  a  priori  condition  that  it  be  either 
of  the  north-central  U.S.  or  the  Gulf  Coast  (U.S.)  types  of 
sediment-covered  stable  regions:  the  first  of  these  regions  in  the 
U.S.  has  been  found  to  have  an  upper  mantle  with  shield  character¬ 
istics  while  the  second  has  an  upper  mantle  with  a  strong  low  velocity 
channel,  and  is  the  prototype  younger  stable  region. 

The  next  three  regions  are  characterized  by  their  high 
tectonic  activity.  We  have  chosen  to  identify  the  mountainous 
collision  zone  between  the  Asian  and  Indian  plates  as  a  single 
province;  this  is  undoubtedly  incorrect  in  detail,  but  is  valid 
enough  in  view  of  our  inability  to  provide  detailed  resolution  of 
smaller  regions  by  our  use  of  long  paths. 

It  should  be  noted  that  an  increase  in  the  density  of  path 
samples  does  not  necessarily  improve  the  resolution  of  structure 
of  a  small  geographic  area.  Consider  the  limiting  case  of  a 
region  which  is  vanishingly  small  in  area.  An  increase  in  the 
number  of  paths  which  transect  this  region  does  not  improve  our 
ability  to  estimate  the  cross-section  under  it,  because  the 
travel-time  spent  by  the  ray  in  this  region  is  a  vanishingly 
small  part  of  the  total  delay.  An  increase  in  the  number  of  paths 
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in  a  small  area  only  increases  the  redundancy  of  the  data.  For  this 
reason  we  had  some  hesitation  regarding  the  possibility  of  determin¬ 
ing  the  structure  under  a  region  as  "small"  as  the  Tibetan  plateau. 
The  Tibetan  plateau  is  indeed  small  in  comparison  with  the  dimen¬ 
sions  of  the  Eurasian  continent,  and  herein  lies  one  of  the  dif¬ 
ficulties  with  the  method  of  very  long  single-station  transverses 
across  the  entire  span  of  the  continent.  But  the  Tibetan  plateau 
contains  one  of  the  important  mysteries  of  modern  Plate  Tectonics, 
namely  the  reason  for  its  great  elevation  and  the  nature  of  the 
architectural  underpinning  that  holds  it  up,  so  we  resolved  to  try 
to  determine  its  structure , fully  anticipating  that  the  error  bias 
in  the  determinations  of  structural  parameters  might  be  large;  it 
became  our  region  4.  To  improve  resolution  here,  we  made  use  of 
shorter  paths  from  two  Hindu  Kush  and  one  Tashkent  earthquakes  plus 
a  Lop  Nor  explosion  all  recorded  at  nearby  stations  to  increase  the 
fraction  of  the  travel-times  spent  in  the  Tibetan  region.  It  can  be 
seen  from  the  table  of  period  ranges  of  the  observations  (Table  3) 
that  the  inability  to  obtain  long-period  information  from  the 
observations  of  the  Tashkent  earthquake  and  the  Lop  Nor  nuclear  ex¬ 
plosion  also  limits  our  ability  to  resolve  deeper  structure  under 
region  4.  The  complete  sampling  of  all  regions  by  the  path-phase 
delays  from  the  set  of  earthquakes  and  explosions  is  shown  in  Fig.  26. 

The  Sinkiang-Mongolia  seismic  zone  is  easily  identified  as  one 
of  the  provinces  in  our  regionalization.  We  have  chosen  to  identify 
the  relatively  stable  aseismic  block  of  Southeast  China 
as  an  additional  province.  The  Hsingtai  earthquakes  of  1968 
occurred  on  the  boundary  between  these  two  regions.  We  have 
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identified  Southeast  China  as  a  separate  region  without  any  a  priori 
assumption  that  it  has  either  a  cross-section  similar  to  those  for 
the  ancient  shields  or  the  younger  stable  continents.  We  prefer  to 
allow  the  inversion  to  yield  this  determination.  Should  the  in¬ 
version  indicate  that  this  region,  or  region  2  for  that  matter,  is 
similar  in  cross-section  to  that  for  some  other  region,  we  may  make 
the  assumption  that  these  regions  are  the  same  and  use  this  information 
to  perform  a  further  simplified  inversion.  Although,  within  them¬ 
selves,  these  regions  encompass  widely  differing  geologic  structures 
and  widely  varying  seismic  activity,  we  have  assumed  that  each  of 
these  regions  is  homogeneous  in  order  to  limit  the  number  of 
parameters  in  the  inversion. 

Three  small  regions  are  treated  specially.  For  the  South 
China  Sea,  we  have  assumed  the  structure  to  be  known,  and  to  be 
that  for  typical  marginal  seas.  The  phase  delays  for  this  region 
are  taken  from  two  single-station  phase  velocity  observations  in 
marginal  seas  obtained  by  Leeds  <1973)  and  from  two  phase  velocity 
determinations  by  the  two-station  method  across  the  Philippine  Sea. 
Phase  velocity  corrections  for  the  Sea  of  Okhotsk  were  obtained 
theoretically  from  a  crustal  structure  given  by  Kosminskaya  et  al 
(1969)  and  derived  from  explosion  work,  in  which  the  Okhotsk 
depression  has  a  25  km  crustal  thickness;  we  have  used  an  oceanic 
mantle  below  this  crust;  These  phase  corrections  have  been  applied 
to  travel  times  for  those  paths  that  traverse  these  two  regions. 
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Finally,  a  small  fraction  of  the  paths  cross  the  Yakutsk-Magadan 
region  of  northeastern  Siberia.  The  maps  show  this  mountainous 
region  to  be  tectonically  different  in  a  significant  way  from 
the  shield  immediately  to  the  west  of  it.  Since  our  samples  of 
this  region  are  all  from  the  Kamchatka  earthquake,  and  are  all 
small  in  length,  we  have  decided,  without  justification,  to 
couple  this  region  to  our  other  mountainous  province,  region  3. 

We  do  not  expect  significant  differences  from  a  region  3  structure 
for  this  part  of  Eastern  Eurasia  to  produce  major  changes  in  the 
inversions  since  the  total  fraction  of  path  length  in  this  region 
is  small. 

The  total  path  length  in  each  region  summed  over  all  event- 
paths  in  given  in  Table  4.  Since  the  phase  shifts  have  not  been 
obtained  over  a  uniform  band  of  periods  for  all  paths,  in  the  last 
column  of  Table  4  we  have  also  reported  an  estimate  of  the 
number  of  samples  in  each  region  by  giving  the  product  of  sample 
path  length  by  number  of  period  estimates  of  phase  shift.  By 
either  method  of  estimation,  the  very  low  fraction  of  sampling 
in  regions  4  and  6  lead  us  to  expect  that  large  uncertainties  in 
the  structure  will  be  obtained  from  the  inversion  for  these 
regions . 

We  have  inverted  the  data  under  the  assumption  that  a  simple 
ray  theory  for  surface  waves  applies,  that  is,  the  phase  shift  for 
a  surface  wave  passing  through  a  given  region  is  computed  as  though 
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TABLE  4 


Reqion 

Path 

Length (km) 

Percentage  of  Total 
Path  Length 

Weighted  Percentage 
Total  Path  Lenqth 

1 

20970 

12.1 

13.0 

2 

53150 

30.6 

30.4 

3 

38411 

22.1 

22.5 

4 

15292 

8.8 

6.9 

5 

37994 

20.7 

20.7 

6 

9979 

5.7 

6.5 

TOTAL 

173796 

100.00 

100.0 

t 
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the  region  were  laterally  infinite  in  extent  and  uninfluenced  tv 
the  presence  of  neighboring  regions,  no  matter  how  close  the  great 
circle  path  approaches  the  regional  boundaries.  The  assumption 
that  diffraction  effects  are  unimportant  is  evidently  untenable 
but  provides  a  basis  for  starting  an  inversion;  we  have  tested 
this  assumption  in  one  of  the  inversions  below. 

The  inversion  proceeds  using  the  method  described  by  Leeds 
et  al .  (1974).  We  calculate  the  phase  travel  time  for  the  n*"^1 

path  and  the  pth  period  as 


6^ 

t  V  *.  s. 

np  =  ^  in  ip 


where  is  the  path  length  of  the  nth  path  in  the  ith  region 

(i=l,...,6)  and  s^  is  the  (theoretical)  phase  slowness  for  the  ith 
region  at  the  pth  period.  The  phase  slownesses  s^  are  functions 
of  the  model  parameters  in  each  region. 


V 

n 


v  (t  -  t  )2 
“  np0  np 


where  t  is  the  observed  travel-time  for  the  nth  path  at  the  pth 

npo 

period.  The  minimization  takes  place  with  respect  to  the  choice 


of  model  parameters. 
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We  obtain  larger  and  larger  variances  in  the  model  para¬ 
meters  as  the  number  of  model  parameters  increases.  We  would 
like  to  be  able  to  solve  for  the  properties  of  the  crust  in  each 
of  the  six  regions.  However,  this  a)  requires  much  precise  data 
at  periods  shorter  than  30  sec.,  b)  increases  the  number  of  model 
parameters  significantly,  and  c)  pushes  our  postulate  of  lateral 
homogeneity  in  each  of  the  regions  to  an  untenable  extreme.  We 
have  therefore  used  model  crusts  for  each  of  the  regions  which  a) 
seem  to  be  plausible  when  compared  with  results  for  similar  parts 
of  the  earth  where  observations  exist  (such  as  locations  typical  of 
regions  1,  2  and  3) ,  and  b)  agree  with  Soviet  refraction  results 
where  available  (Kosminskaya,  et  al,  1969;  Sollogub,  1969).  When 
large  residuals  were  encountered  at  short  periods,  such  as  in  the 
case  of  region  3  (Alpide-Himalayan  belt)  and  region  4  (Tibetan  plateau) , 
we  were  obliged  to  introduce  more  low-velocity  material  into  the  crust. 
This  was  done  by  keeping  crustal  velocities  fixed  and  increasing 
crustal  thickness.  In  these  two  cases,  this  procedure  leads  to 
extraordinarily  thick  crusts.  It  should  be  realized  that  these 
model  crustal  thicknesses  are  consequences  of  the  procedures 
used;  if  we  had  chosen  to  lower  the  crustal  velocities,  the 
thicknesses  would  have  been  less.  We  have  used  a  crustal  thick¬ 
ness  of  70  km  in  the  Tibetan  Plateau  (region  4)  a  value  not  in¬ 
consistent  with  other  recent  estimates  (Chun  and  Yoshii,  1977). 

The  crustal  models  we  have  used  are  listed  in  Table  5  and  are 
assumed  to  be  fixed  in  the  inversions. 


IABlL  i, 
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Model  parameters  for  inversion  of  phase  velocity  data 


CRUST 

LID 

CHANNEL 

SUBCHANNEL 


Thickness  Depth  B(km/sec)  a(km/sec) 

-  0 


(different  crustal  models  for  each  region) 

-  (fixed) 

VAR  4.65  8.17 


VAR  7.80 


VAR  4.80  8.80 

-  420( fixed) 


94 


94 


63 


138 


139 


104 


514 

608 

691 

809 

948 

1052 


5.128 

5.283 

5.419 

6.172 

6.266 

6.351 


9.609 

9.781 

5.902 

10.552 

11.238 

11.392 


l 


L 


a 


p(gm/cm3) 

3.45 

3.45 

3.65 

3.806 

3.934 

4.051 

4.417 

4.505 

4.579 


TABLE  5  (cont'd) 
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We  have  parameterized  the  mantle  into  a  lid,  channel  and 
subchannel;  each  of  these  layers  in  a  given  region  is  taken  to  be  homo¬ 
geneous.  The  subchannel  region  terminates  at  a  depth  of  420  km. 

Below  this  depth,  we  place  a  standard  lower  mantle  platform  under 
all  regions,  indicated  in  Table  5.  The  unknowns  in  the  inversions 
will  be  the  shear  velocities  in  each  of  the  three  upper  mantle 
layers  and  the  depths  of  the  interfaces  at  top  and  bottom  of  the 
second  ("channel")  layer.  The  number  of  unknowns  is  thus  30, 
five  in  each  of  the  six  regions.  This  figure  far  exceeds  the  number 
of  degrees  of  freedom  in  the  data.  Thus  some  reduction  in  the 
number  of  unknowns  has  to  be  made,  by  some  assumptions  which  are 
geophysical  in  character. 

Inversion  1 

In  a  first  attempt  at  inversion,  we  have  fixed  the  lid  S-wave 
velocity  at  4.65  km/sec  and  the  subchannel  S-wave  velocity  at  4.8 
km/sec  in  order  to  reduce  the  number  of  degrees  of  freedom.  The 
value  of  4.65  km/sec  for  the  lid  arises  frequently  in  inversions 
for  other  parts  of  the  world.  For  one  case  in  which  a  4.65  km/sec 
lid  was  not  observed,  namely  for  the  Western  United  States  (Biswas 
and  Knopoff,  1974)  in  which  the  subcrustal  material  has  S-wave 
velocity  around  4.3  km/sec,  we  were  able  to  assume  that  a  model  with  a 
4.65  km/sec  lid  with  zero  lid  thickness  was  accepted  by  the  inversion 
and  that  the  4.3  km/sec  value  represents  channel  material  rising 
almost  to  the  base  of  the  crust.  Should  the  lid  velocity  in  some 
region  be  less  than  4.65  km/sec  in  reality,  then  crustal  thick¬ 
nesses  can  be  reduced.  A  similar  comment  can  be  made  about  the 
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subchannel  velocity:  if  the  subchannel  velocity  should  turn  out  to 
be  less  than  4.8  km/sec  in  some  regions,  channel  thicknesses  will 
be  reduced. 

This  parameterization  thus  includes  only  two  adjustable  model 
parameters  for  each  region.  These  are  the  lid  thickness  and  the 
channel  S-wave  velocity.  A  thirteenth  parameter  is  the  sub¬ 
channel  thickness,  which  is  presumed  to  be  uniform  across  the 
entire  continent  and  hence  has  the  same  value  under  each  region. 
Since  the  crustal  thickness  and  the  depth  to  the  420  km  interface 
are  fixed,  the  parameterization  of  lid  and  subchannel  thicknesses 
is  equivalent  to  a  parameterization  of  the  depth  below  the  sur¬ 
face  of  the  top  and  bottom  of  the  channel.  This  parameterization 
has  13  degrees  of  freedom. 

After  adjustment  of  the  crust  by  the  procedures  described 
above  (with  interpretation  of  crustal  parameters  according  to 
the  remarks  above) ,  the  parameterization  and  cross  sections  used 
in  a  linearized  inversion  procedure  are  given  in  Table  5.  The 
superficial  sedimentary  layer  that  is  introduced  in  the  crusts 
of  regions  3  and  4  is  designed  to  reduce  the  residuals  at  the 
shortest  periods. 
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The 

starting  values 

of  the 

thirteen 

parameters  in 

this  linearized 

inverse 

are (in  the  usual 

units) 

: 

Table  6 

Region 

PCH 

hLID 

hSUB 

Depth  to 
top  of 
channel 

Depth  to 
bottom  of 
channel 

1 

4.51 

110 

;■ 

155 

", 

2 

4.39 

113 

’ 

158 

3 

4 . 30 

54 

150 

119 

270 

4 

4.29 

90 

167 

5 

4.08 

92 

i 

137 

6 

4.38 

103 

148 

The  thirteen  parameters  in  the  rectangular  box  in  Table  o  are  those  var¬ 
ied  in  the  inversion.  The  standard  errors  of  the  travel-times  in  the 
inversion  were  taken  to  be  the  same  as  those  used  by  Leeds  et  al . 

(1974),  namely  o  =  Max  (7,  0.1T)  sec.  These  estimates  were  later 
found  to  be  too  large  and  required  modification 

In  the  inversion,  an  iteration  process  has  been  used  in  which 
the  matrix  of  partial  derivatives  G  was  recalculated  whenever  we 
moved  into  a  new  portion  of  parameter  space.  This  process  was  con¬ 
tinued  until  we  obtained  convergence  of  the  variable  model  parameters. 

The  thirteen  eigenvalues  of  the  product  matrix  of  partial  derivatives 
T 

G  G,  which  were  obtained  in  the  final  stage  of  the  iteration  process, 
are  given  in  Table  7.  Each  eigenvalue  corresponds  to  an  eigenvector 
which, in  every  case  (except  nos.  5  and  6), points  in  a  direction  close 
to  one  of  the  thirteen  parametric  axes.  Thus,  each  eigenvector  can 
be  said  to  be  a  discriminant  for  each  of  the  thirteen  degrees  of 
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TABLE  7 

Eigenvalue 


1. 

13.68 

2. 

8.43 

3. 

6.85 

4. 

3.56 

5. 

2.87 

6. 

2.59 

7. 

1.84 

00 

1.61 

9. 

1.14 

10. 

1.04 

11. 

0.56 

12. 

0.36 

13. 

0.21 

Model  parameter  most  closely  resolved 
by  eigenvector 

~CH2 

6CH5 

BCH3 

eCHl 

SCH6  &  hSUB 

feCH6  &  hSUB 

6CH4 

hLID5 

hLID2 

hLID3 

hLID6 

hLIDl 

hLID4 
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freedom  in  the  model.  The  eigenvectors  corresponding  to  the 
largest  eigenvalues  give  reliable  structural  information  con¬ 
cerning  the  model  parameters  closest  to  them,  and  the  eigen¬ 
vectors  corresponding  to  the  smallest  eigenvalues  give  little 
information  about  their  corresponding  parameters. 
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From  Table  7  it  can  be  seen  that,  in  general,  the  data  yield 
the  greatest  information  about  the  channel  velocities  and  the 
least  information  about  the  lid  thicknesses.  In  general,  the 
data  contain  more  information  about  regions  2,  3,  and  5  than 
the  the  other  three  regions.  This  is  clearly  correlated  with 
the  higher  percentages  of  total  path  length  which  sample  regions 
2,  3,  and  5.  (For  complete  details,  see  Table  4  and  Fig.  26). 

The  variance  of  each  model  parameter  decreases  as  successive 
models  approach  the  fine  one,  which  demonstrates  the  convergence 
of  our  inversion  procedure.  Fig.  27  gives  our  results  for  this 
inversion:  the  "best"  model  and  the  corresponding  standard 

deviations  of  each  model  parameter.  We  see  that  the  upper  mantle 
structure  for  regions  l  and  2  are  very  similar  although  the  lid 
in  region  1  is  somewhat  thicker.  The  channel  shear-wave  velocity 
for  region  1  is  slightly  less  than  that  in  region  2  ,  but  the  un¬ 
certainty  of  the  model  parameter  in  region  1  is  rather  large  and 
in  any  case,  the  velocity  contrast  to  the  lid  in  these  regions 
is  rather  smaller  than  one  would  like  to  assert  that  a)  a  channel 
is  indeed  present  and  b)  that  a  channel  with  velocities  low 
enough  to  require  partial  melting  is  present.  Thus,  the  existence 
of  a  low-velocity  channel  in  region  1  is  uncertain.  The  most  strik¬ 
ing  result  of  this  inversion  is  the  very  thin  lid,  and  moderate 
shear  wave  velocity,  in  the  channel  for  regions  3  and  5;  as 
in  the  case  of  regions  1  and  2,  regions  3  and  5  have  strikingly 
similar  upper  mantles.  These  latter  two  regions  are  those  which 
are  tectonically  active  and  have  high  seismicity.  Region  4 
is  that  in  which  the  Eurasian  and  Indian  continents  are 
in  collision.  Although  the  upper  mantle  structure  for  this 
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region  appears  to  be  very  similar  to  those  of  regions  1  and  2, 
the  uncertainties  in  the  model  parameters  are  large.  This  is  due 

1)  to  the  low  percentage  of  total  path  length  in  this  region,  and 

2)  to  the  fact  that  most  of  the  paths  which  sample  this  region 
have  only  short-period  information.  Region  6  has  a  thick  lid  and 
a  pronounced  low-velocity  channel,  but  the  standard  deviations  are 
rather  large;  additional  assumptions  may  help  to  improve  the 
resolution  of  the  structural  parameters  in  this  region. 

Inversion  2 

In  a  second  attempt  at  inversion  we  have  made  the  assumption 
that  the  chemical  composition  of  the  lids  is  everywhere  the 
same  across  Eurasia  (as  before)and  that,  in  this  case,  the  chemical 
composition  of  the  channels  is  similarly  the  same  across  Eurasia. 

We  have  allowed  the  lid  velocity  to  be  adjusted  in  this  case;  it 
was  fixed  in  the  preceding  case.  In  this  case  the  number  of  para¬ 
meters  in  the  fit  is  nine;  the  depth  to  the  lid-channel  interface 
in  each  of  the  six  regions,  the  depth  to  the  channel  -  subchannel 
interface,  and  the  S-wave  velocities  in  the  lid  and  channel; 
the  latter  three  parameters  are  are  constants  across  the 
entire  region.  The  errors  o  are  assumed  to  be  the  same  as 
before . 

The  results  of  this  inversion  are  shown  schematically  in  Fig.  28. 
The  result  for  lid  velocity  was  4.56+. 01  km/sec  and  for  channel 
velocity  4.34+. 02  km/sec.  The  results  for  lid  thicknesses  show  striking 
similarity  among  regions  1  ,  2,  4  on  the  one  hand  and  3,  5,  and  6  on  the  other. 


FIG.  26.  Propagation  path  in  each  subregion  of  the 

Eurasian  continent.  The  percentage  of  total 
path  length  in  each  region  is  shown  in  Table  2. 

FIG.  27.  Final  result  for  linearized  inversion  of  13 
variable  parameters  in  the  model  and  their 
corresponding  standard  deviations. 

FIG.  28.  Schematic  upper  mantle  cross-sections  obtained 


in  inversion  2. 
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This  observation  lends  support  to  the  interpretation  of  the  first 
inversion,  that  che  sediment-covered  platforms  of  Russia  and  Siberia 
are  largely  sediment-covered  shields.  The  result  for  Tibet  we 
have  obtained  is  that  the  upper  mantle,  to  great  depth,  is 
largely  shield-like  in  character.  We  interpret  this  to  mean, 
if  it  should  be  substantiated,  that  the  Indian  Shield  has  been 
emplaced  beneath  the  Tibetan  crust  during  the  collision  of  the 
plates  and  that  the  mantle  beneath  Tibet  is  presently  at 
temperatures  well  below  the  melting  point  to  great  depth,  i.e. 
the  mantle  under  the  Tibetan  plateau  is  relatively  cool  com¬ 
pared  with  tectonically  active  collision  zones  such  as  region  3. 

On  the  other  hand  regions  3,  5,  and  6  all  exhibit  the  develop¬ 
ment  of  a  marked  low-velocity  channel  at  shallow  depths  in  the 
mantle,  which  implies  that  a  significant  zone  of  partial  melting, 
perhaps  200  km  thick,  is  present.  All  these  three  regions  are 
thus  presumed  to  be  tectonically  active,  despite  the  low  seismicity 
of  region  6  (the  large  recent  earthquakes  in  China  occurred 
in  the  northern  part  of  region  6).  The  presence  of  a  small  channel 
at  great  depth  is  not  considered  to  be  significant:  this  may  be 
an  artifact  of  the  inversion,  due  to  improper  choice  of  parameters 
at  shallower  depths.  These  deep  channels  are  probably  absent 
but  we  cannot  be  absolutely  certain. 

Inversion  2.1 

We  have  made  a  variation  on  inversion  2  by  enlarging  the 
extent  of  the  Baltic  shield  in  the  preceding  inversion,  as  shown 
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in  Fig.  29  according  to  a  hint  in  the  tectonic  map  of 
Khain  and  Muratov  (1969).  As  might  have  been  expected,  the 
results  of  the  preceding  inversion  are  not  significantly  changed 
in  view  of  the  result  of  inversion  2,  namely  that  regions  1  &  2 
are  virtually  identical  in  cross-section.  Hence  enlarging 
region  1  at  the  expense  of  region  2  cannot  be  expected  to  pro¬ 
duce  a  significant  change  in  the  results. 

A  tabulation  of  the  results  of  this  inversion  is  as  follows, 
and  is  listed  only  for  purposes  of  comparison  with  the  other 
inversions : 


^LID 

= 

4.56+0.01  km/sec 

*CH 

= 

4.34+0.02  km/sec 

hLID 

1' 

= 

206+25  km 

hLID 

2' 

= 

229+29 

hLID 

3 

= 

99+9 

hLID 

4 

= 

246+56 

hLID 

5 

= 

79  +  7 

hLID 

6 

= 

69  +  7 

hCH-SUB 

= 

274  +  8 

Inversion  2.2 

In  each  of  the  two  preceding  inversions  we  have  obtained 
the  result  that  the  upper  mantle  of  both  regions  1  and  2  are 
remarkably  similar.  In  both  cases  we  have  a  deep  continental 
root  to  essentially  the  same  depths.  We  have  therefore  made 
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29.  Modification  of  regionalization  of  Fig.  25 
by  enlargement  of  region  1  in  Baltic  shield 
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the  assumption  that  these  two  regions  are  indistinguishable 
in  a  new  inversion  2.2.  We  call  the  coalescence  of  regions 
1  and  2,  a  new  region  1A  (Fig.  30);  there  is  now  no  separate 
region  2. 

As  might  be  expected  the  results  of  the  inversions  are 
not  significantly  changed  from  the  earlier  explorations.  Indeed 
the  results  for  inversion  2.2  are  identical  in  all  respects 
with  that  of  inversion  2.1 


6lid 

= 

4.56+. 01  km/sec 

6ch 

= 

4.34+. 02  km/sec 

hCH-SUB 

= 

274+8  km 

hLID(IA) 

= 

217+18  km 

hLID(3) 

= 

99+9 

hLID(4) 

= 

246+55 

hLID(5) 

= 

79+9 

hLID  (6) 

69+7 

all  depths  h  are  measured  from  the  surface. 

Inversion  3 

A  rather  annoying  aspect  of  the  data  concerns  the  fact  that 
the  ph^se  travel-time  residuals  from  any  of  the  preceding  models 
are  not  normally  distributed.  There  is  an  unacceptably  large 
number  of  residuals  between  2 a  and  4c;  this  result  has  been 
verified  by  a  x2  test.  We  have  deleted  several  of  the  data  with 
large  residuals  and  have  proposed  a  new  data  set  with  a  value  of 
X2  which  places  the  new  data  set  within  acceptable  limits  of  a 
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FIG.  30.  Modification  of  regionalization  of  Fig.  25 


by  fusion  of  regions  I  and  2. 
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a  normal  distribution  with  a  new  value  of  a  having  2/3  the 
former  value: 

a  =  2/3  Max  (7  sec,  0.1T) 

We  have  repeated  inversion  2.2  with  the  revised  data 
set  and  the  new  estimate  of  errors  and  have  obtained  the  follow¬ 
ing  results: 


PLID 

4.57+. 01  km/sec 

bch 

= 

4 . 35+ . 01  km/sec 

hCH-SUB 

= 

276+5  km 

hLID(IA) 

= 

204+12 

hLID(3) 

= 

99  +  7 

hLID(4) 

= 

256+41 

hLID (5) 

= 

59  +  4 

hLID(6) 

= 

55  +  5 

The  shear  velocities  in  the  lid  and  channel  are  unchanged  from  the 
preceding  inversion  2.2.  But  there  have  been  some  significant 
changes  in  the  lid  thickness  of  some  of  the  regions.  Regions  1A 
and  4  continue  to  have  upper  mantles  consistent  with  ancient  (cold) 
shield  models.  However,  regions  5  and  6  now  have  very  thin  lids, 
implying  high  heat  flow  and  the  presence  of  strong  tectonic  pro¬ 
cesses  . 

Inversion  3.1 

To  test  the  validity  of  the  linearized  inversion  procedures 
especially  in  view  of  our  observations  that  the  final  results  are 
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strongly  dependent  on  smoothing  procedures  in  the  model  space, 
we  have  applied  the  power  of  our  non-linear  Hedgehog  program 
whose  results  are  independent  of  smoothing  in  the  model  space. 

The  data  set  is  the  reduced  data  set  of  Inversion  3  with  a 
normal  distribution  of  residuals  relative  to  model  2.2.  We 
have  further  fixed  the  channel-subchannel  interface  at  the  value 
given  by  inversion  3.  The  remaining  seven  parameters  were  ex¬ 
plored  in  the  space  given  by  Table  6. 


TABLE  6 

Parameters  and  Range  of  Search  in  Hedgehog  Inversion 
for  Inversion  3.1 


Starting  Value 

Step  Size 

Lower  Limit 

Upper  Limit 

PI 

4.57 

0.1 

4.47 

4.67 

P2 

4.35 

0.1 

4.25 

4.45 

P3 

204 

30 

144 

264 

P4 

99 

30 

69 

189 

P5 

256 

30 

136 

256 

P6 

59 

30 

59 

179 

P7 

55 

30 

55 

175 

PI: 

Lid  shear  wave 

velocity 

in 

Km/Sec 

P2 : 

Channel  shear 

wave  velocity  in  Km/Sec 

P3 : 

Lid  thickness 

of  region 

1A 

in  Km 

P4 : 

Lid  thickness 

of  region 

3 

in  Km 

P5 : 

Lid  thickness 

of  region 

4 

in  Km 

P6 : 

Lid  thickness 

of  region 

5 

in  Km 

P7 : 

Lid  thickness 

of  region 

6 

in  Km 
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The  results  of  the  inversion  gave  the  67  acceptable  solutions, 
acceptable  within  1.4a  under  the  postulate  that  there  is  equal 
tradeoff  between  the  effects  of  uncertainties  in  the  model  space 
and  fit  to  the  data  (Table  7).  The  first,  and  most  obvious 
result  is  that  the  lid  velocity  was  accepted  to  be  4.57  km/sec; 
no  other  values  are  acceptable.  The  channel  velocities  accepted 
were  only  4.25  and  4.35  km/sec.  Although  the  lid  thicknesses 
for  the  shields  of  region  1A  could  vary  between  174  and  234  km 
and  for  the  Tibetan  plateau  (region  4)  between  226  and  256  km 
a)  under  some  circumstances  both  lids  could  have  the  same  thick¬ 
ness  (models  14,  30,  31,  etc.)  and  b)  under  no  circumstances 
could  the  Tibetan  plateau  have  a  lid  which  was  as  thin  as  (say) 

30  km,  a  result  which  would  have  implied  a  mantle  appropriate 
to  a  tectonically  active  region,  i.e.,  one  with  a  well-developed 
high-contrast  low  velocity  channel.  Finally,  although  the 
linear  inverse  of  inversion  3  implied  a  difference  in  lid 
thickness  between  region  3  on  the  one  hand  and  regions  6  and  7 
on  the  other,  the  non-linear  model-independent  inverse  gave  some 
solutions  in  which  the  lid  thicknesses  are  such  that  the  lid 
channel  interface  is  roughly  at  a  common  depth  below  the 
surfaces  of  all  three  regions  (e.g.  solutions  35,  55,  23,  37, 

46,  57).  We  conclude  from  these  inversion  results  that  cold 
shield  properties  extend  to  great  depth  under  regions  1A  and  4 
while  regions  3,  5,  6  may  have  similar  properties  (at  least 
we  cannot  discount  this  result)  with  upper  mantles  characteristic 
of  young  active  regions.  What  is  remarkable  is  that  region  6 
is  coupled  together  with  the  more  obvious  active  regions. 


Table  7 
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Successful  solutions  in  Hedgehog  inversion  3.1 


t'l 

i’2 

Pi 

P4 

i'5 

16 

F7 

1 

4.57 

4.35 

204 

99 

256 

59 

55 

2 

4.57 

4.35 

174 

99 

256 

5 

55 

3 

4.57 

4.35 

234 

99 

256 

59 

55 

4 

4.57 

4.35 

204. 

69 

256 

59 

55 

5 

4.57 

4.35 

204 

129 

256 

59 

55 

6 

4.57 

4.35 

204 

GO 

226 

59 

55 

7 

4. 57 

4.35 

204 

99 

256 

39 

55 

0 

A.  57 

4.  30 

204 

00 
✓  / 

256 

59 

65 

9 

4.57 

iu55 

174 

69 

256 

59 

55 

10 

4.57 

4.35 

174 

129 

256 

59 

55 

1 1 

4.57 

4. 35 

234 

69 

256 

59 

55 

12 

4.57 

4.35 

234 

129 

256 

59 

55 

13 

4.57 

4.35 

174 

99 

226 

59 

55 

14 

4.57 

4.35 

234 

go 

226 

59 

55 

15 

4.57 

4.35 

174 

99 

256 

89 

55 

16 

4.57 

4.35 

234 

99 

256 

89 

55 

17 

4.57 

4.35 

174 

99 

256 

59 

85 

IS 

if.  57 

4.35 

234 

99 

256 

59 

85 

19 

4.57 

4.35 

204 

69 

226 

59 

55 

20 

4.57 

4.35 

204 

129 

226 

59 

55 

21 

4.57 

4. 35 

204 

69 

256 

89 

55 

22 

4.57 

4.35 

204 

129 

256 

69 

55 

23 

4.57 

4. 35 

204 

69 

256 

59 

85 

24 

4.57 

4. 35 

204 

129 

256 

59 

35 

25 

4.57 

4.35 

204 

99 

226 

89 

55 

26 

4.57 

4.35 

204 

00 

X1  y 

226 

59 

85 

27 

4.57 

4.35 

204 

99 

256 

89 

35 

28 

4.57 

4.35 

174 

69 

226 

59 

55 

29 

4.57 

4.35 

174 

129 

226 

59 

55 

50 

4.57 

4.35 

234 

69 

226 

59 

55 

31 

4.57 

4.35 

234 

129 

226 

59 

55 

32 

4.57 

4.35 

174 

69 

256 

89 

55 

33 

4.57 

4.35 

174 

129 

256 

89 

55 

34 

4.57 

4.35 

234 

69 

256 

89 

55 

L 
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i'1 

r  2 

r  J1' 

4.57 

4.35 

17'. 

36 

4.57 

4.35 

174 

37 

4. 57 

4.35 

234 

36 

4. 57 

4.35 

234 

39 

4.5? 

4.35 

174 

40 

4.57 

4.35 

234 

41 

4.5? 

4.35 

174 

42 

4. 57 

4. 35 

234 

43 

4.57 

4.35 

174 

44 

4.57 

4*35 

204 

45 

4.57 

4.35 

204 

46 

4.5? 

4.35 

204 

47 

4. 5? 

4.35 

204 

43 

4. 57 

4.35 

204 

49 

4.57 

4. 35 

204 

50 

4.57 

4. 25 

234 

51 

4.57 

4.25 

204 

52 

4.57 

4.35 

174 

53 

4.57 

4.35 

174 

54 

4.57 

4.35 

234 

55 

4.57 

4.35 

174 

56 

4.5? 

4.35 

174 

57 

4.5? 

4.35 

234 

56 

4.5? 

4.  J^5 

234 

59 

4.57 

4.35 

174 

60 

4.57 

4.35 

234 

61 

4.57 

4.35 

174 

62 

4.57 

4.35 

234 

63 

4.57 

4.35 

204 

64 

4.57 

4.25 

234 

b5 

4.57 

4.35 

174 

66 

4.57 

4.55 

234 

67 

4.57 

4.25 

234 

24 

i'5 

i-6 

i7 

69 

256 

59 

o5 

129 

256 

59 

85 

69 

256 

59 

85 

129 

256 

RQ 

85 

oo 
✓  ^ 

->  'jjf 

89 

55 

99 

226 

c° 

55 

99 

226 

59 

85 

99 

226 

59 

85 

99 

256 

69 

65 

69 

226 

89 

55 

129 

226 

89 

55 

69 

226 

59 

85 

129 

226 

59 

85 

69 

256 

89 

35 

99 

226 

89 

65 

129 

256 

69 

55 

129 

256 

89 

85 

69 

226 

89 

55 

129 

226 

89 

55 

69 

226 

89 

55 

69 

226 

59 

65 

129 

226 

59 

85 

69 

226 

59 

65 

129 

226 

59 

85 

69 

256 

69 

85 

69 

256 

69 

85 

99 

226 

89 

85 

99 

22  6 

89 

85 

69 

226 

69 

85 

129 

256 

89 

35 

69 

226 

69 

6; 

69 

226 

89 

85 

129 

226 

89 

85 

L 
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Inversion  4 

To  test  the  importance  of  diffraction,  we  have  deleted 
from  the  data  set, phase  travel  times  for  those  paths  which  lie 
close  and  parallel  to  a  contrast  between  two  dissimilar  zones. 
The  value  of  a  continues  to  be  taken  as  in  the  result  of 
Inversion  3.  The  linear  inverse  using  the  same  regionalization 
and  parameterization  as  inversion  3  for  the  results: 


^LID 

= 

4.58+. 01 

km/sec 

6CH 

= 

4 . 36+ . 02 

km/sec 

hLID (1A) 

= 

181+13 

km  (below 

surface) 

hLID(3) 

= 

101+16 

it 

" 

hLID(4) 

= 

204+35 

tt 

ii 

hLID(5) 

= 

54+7 

" 

m 

hLID(6) 

= 

54+7 

II 

ti 

hCH-SUB 

= 

268  +  6 

II 

" 

We  detect  no  significant  changes  from  the  results  of  inversion  3 

Inversion  5 

Finally,  since  the  Tibetan  structure  appears  to  be 
associated  with  the  collision  of  the  Indian-Asian  plates,  we 
have  incorporated  the  Himalayan  part  of  region  3  into  region  4, 
and  re-analyzed  the  linear  inversion  with  the  data  and  errors 
as  in  inversion  3.  A  map  of  the  new  ’-egionalization  is  shown 
in  Figure  31. 
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The  results  of  the  inversion  are 


blid 

— 

4.58+. 01  km/sec 

bch 

= 

4. 34+.01 

hLlD(lA) 

= 

194+10  km 

hLID(3) 

= 

9  7  +  8 

hLID (4) 

= 

207+21 

hLID (5) 

= 

67+4 

hLID (6) 

= 

58  +  5 

hCH-SUB 

= 

27  3+5 

Again  we  observe  no  major  change  in  the  structure. 

We  conclude  from  all  these  tests  that  the  stable  sediment- 
covered  regions  of  the  USSR  are  probably  stable  preCambrian 
shields  covered  by  sediments,  that  the  Tibetan  plateau  is  under¬ 
lain  by  relatively  cold  shield  mantle  material  with  no  major 
low-velocity  zone  at  even  moderate  depths  that  might  be  ex¬ 
pected  of  tectonically  active  zone,  or  of  young  stable  regions, 
that  the  mountainous  collision  zone  between  the  Asian  and 
Indian  plates  have  a  prominent  low  velocity  zone  at  moderate 
depths,  that  the  Sinkiang-Mongolian  active  seismic  zone  has  a 
similar  structure  and  that  South  Eastern  China  is  also  a 
region  of  tectonic  activity  as  indicated  by  the  similarity  of 
its  uppper  mantle  structure  to  the  other  two  seismically 
active  regions. 


46 


Arctic  Study.  We  have  'Measured  fundamental  mode  Rayleigh 
waves  over  a  number  of  paths  crossing  the  Arctic  Ocean.  We 
have  used  as  sources  four  earthquakes  whose  focal  parameters 
a  re  : 


4.  Lena  River 

Aug . 

25,  1964 

13:47:20.6 

78 . 2  °N 

126.6 

°E 

10.  Eastern  Aleutians(l) 

Feb . 

6,  1965 

01  :40  :33.2 

5  3 . 2  °  N 

161.9 

°W 

11-  Eastern  Aleutians{2) 

Feb . 

6,  1965 

16:50:23.6 

53. 3°.N 

161  .8 

°W 

13.  Western  Aleutians 

Feb. 

6,  1965 

04:02:53 

52 . 1 °N 

175.7 

°E 

For  each  of  these  earthquakes  we  have  obtained  initial  phases 
either  from  the  fault  plane  solution  or  from  the  radiation  pat¬ 
tern  for  Rayleigh  waves  (as  in  the  case  of  the  Lena  River  dis¬ 
cussed  above).  We  have  obtained  phase  velocities  by  the  single¬ 
station  method  for  seven  paths  crossing  the  Arctic  over  the 
period  range  50  to  208  sec.  The  paths  are  shown  in  Figure  32  . 

It  can  be  seen  that  none  of  these  are  purely  oceanic  paths. 
The  shaded  area  outlines  out  estimate  of  the  boundary  between 
the  continental  shelf  and  the  deep  ocean  basins.  The  fraction 
of  the  geometrical  path  that  each  event  has  in  the  oceanic  part 
is  as  foil ows  : 


Event 

Total  path  length 

Oceanic  length 

Fraction  oceanic 

1  .  Lena  -  ESK 

4816  km 

1904  km 

.40 

2.  Lena  -  KTG 

3386 

1823 

.54 

3.  East  Aleut(l)  -  KON 

7483 

4603 

.62 

4.  East  Aleut(l)  -  KTG 

5933 

916 

.15 

5.  East  Afeut(2)  -  KEV 

6347 

2334 

.37 

6.  East  Aleut(2)  -  ESK 

7818 

2714 

.35 

7.  West  Aleut.  -  KEV 


6263 


844 


.13 
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FIG.  32.  Map  of  Arctic  region  showing  all  propagation 
paths  used  in  dispersion  and  regionalization 
analysis . 
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To  reduce  the  available  data  to  information  regarding 
purely  oceanic  paths,  we  have  decided  to  use  the  phase  veloci¬ 
ty  data  for  the  Lena  River  event  recorded  at  KEV  (see  discus¬ 
sion  above  for  Eurasia)  as  a  typical  continental  value  and  to 
subtract  these  values  ,  for  the  appropriate  path  length  con¬ 
tribution,  from  the  phase  delays  observed  for  the  above  7 
path-events.  Unfortunately,  the  two  events  East  Aleut.(l)  -  KTG 
and  West  Aleut. -KEV  have  such  small  parts  of  their  total  path 
that  are  oceanic  that  we  are  subtracting  two  numbers  of  com¬ 
parable  size  and  the  result  is  quite  unstable.  The  unreliabil¬ 
ity  of  the  oceanic  phase  delay  results  for  these  two  cases  has  obliged 
us  to  exclude  them  from  our  data  set.  Accordingly,  we  have  in¬ 
vestigated  the  inversions  of  the  phase  velocity  results  for 
the  five  remaining  paths.  The  relevant  data  are  given  in  Table  8: 


TABLE  8 

(Pure)  Oceanic  Phase  Velocities  (km/sec) 


T( sec ) 

Lena -E5K 

Lena -KTG 

E  .  A 1 (1 ) -K0N 

E  .  A 1 ( 2 ) -KEV 

E.A1 ( 2 ) -ESK 

208 

4.64 

(4.83) 

4.62 

4.55 

4.52 

192 

4.52 

(4.61) 

4.51 

4.42 

4.42 

167 

4.38 

4.32 

4.38 

4.22 

4.28 

139 

4.21 

4.14 

4.24 

4.09 

4.15 

119 

4.08 

4.06 

4.12 

4.04 

4.11 

100 

4.01 

3.99 

4.08 

3.96 

4.07 

69 

3.95 

3.90 

4.00 

3.91 

4.00 

50 

3.  91 

3.82 

3.98 

3.89 

3.92 

With  so  few  data,  we  have  not  been  able  to  regionalize  the 
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the  deep  Arctic;  the  number  of  degrees  of  freedom  in  the  data 
is  too  small.  The  best  we  can  do  is  to  consider  the  deep  Arctic 
as  a  single  province  and  investigate  the  consequences  of  invert¬ 
ing  an  "average"  phase  velocity  for  the  region.  The  average 
phase  velocity  is  obtained  from  the  above  table  by  weighting  by 
the  oceanic  path  length  in  each  case.  The  result  is  (omitting 
the  quantities  in  parentheses)  : 


T (sec) 

c (km/sec) 

208 

4.59 

192 

4.47 

167 

4.32 

139 

4.18 

119 

4.10 

100 

4.03 

69 

3.97 

50 

3.92 

These  results  can  be  compared  with  those  obtained  for  Pacific 
paths  by  Leeds  (1973)  from  inversion  of  trans-Pacific  phase 
velocity  data  by  methods  similar  to  those  described  above  for 
trans-Eurasian  paths.  The  pure-age  phase  velocities  for  the 
Pacific  can  be  derived  from  the  cross-sections  resulting  from 
the  inversions;  these  are  shown  in  Fig.  33  for  Pacific  ages 
0-10  my,  20-40  my,  85-110  my.  The  Arctic  data  points  are  shown 
as  circles.  The  Arctic  cross-section  averages  out  to  about  a  30  my 
Pacific  structure.  According  to  the  model  of  Parker  and  Oldenburg 
(1973),  the  lid  thickness  as  a  function  of  the  age  is 
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z(t)  =  9.4t2  km,  with  t  the  age  in  my.  Thus,  assuming  that  the 
lid  and  channel  S-wave  velocities  are  of  the  same  order  as  in 
the  Pacific,  the  age  at  which  the  Arctic  began  to  open  is  cal¬ 
culated  to  be  about  70  my  (before  present) ,  an  unexpectedly  small 
quantity. 

Theoretical  seismograms.  The  main  thrust  of  our  work  with 
time  series  synthesis  has  been  directed  toward  improving  the 
efficiency  of  existing  computational  techniques.  This  improve¬ 
ment  is  required  to  permit  us  to  extend  the  information  contained 
on  the  theoretical  seismograms  down  through  the  period  range 
covered  by  the  long-period  instruments  of  the  WWSSN.  Although  up 
to  the  present  time  we  have  concentrated  on  laterally  homogeneous 
structures,  in  all  other  respects  our  models  of  the  earth  have  been 
highly  realistic:  approximately  200  layers  are  being  used  to  model 
the  radial  heterogeneity  of  the  crust-mantle  system  of  a  spherical 
earth,  and  the  intrinsic  attenuation  is  included. 

A  summary  of  the  general  methods  we  have  applied  in  our 
computations  is  given  by  Kausel  and  Schwab  (1973),  and  Knopoff 
et  al  (1974).  An  elaboration  of,  and  certain  justifications 
for  these  procedures  have  recently  been  given  by  Schwab  and 
Kausel  (1976) .  A  recent  contribution  by  Calcagnile  et  al  (1976) , 
also  completed  under  this  contract,  is  also  pertinent  here  when 
only  the  surface-wave  portion  of  the  thoeretical  seismogram  is 
desired. 

The  initial  development  of  the  algorithm  and  certain  pro¬ 
gramming  improvements,  which  were  carried  out  under  the  present 
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FIG.  33.  Results  of  phase  velocity  dispersion  for  Arctic 
region  (circles).  "Standard"  phase  velocity 
curves  for  Pacific  age  strips  are  shown  for 
comparison  (Leeds,  1973).  Because  the  Pacific 
spreads  at  a  different  rate  than  the  Arctic, 
this  gives  a  different  age  for  the  Arctic  (see 
text) . 
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contract,  are  contained  in  Nakanishi,  Schwab  and  Kausel  (1976), 
and  Nakanishi,  Schwab  and  Knopoff  (1976).  These  manuscripts 
formed  the  necessary  bridge  between  first  generation  and  second 
generation  dispersion  programs.  First  generation  techniques 
are  based  on  full,  detailed  specification  —  a  structure,  a 
specific  mode,  and  a  specific  period  —  from  which  a  single 
phase  velocity  at  a  time  is  sought;  second  generation  dispersion 
computations  begin  with  only  the  structure  and  the  mode  specified, 
and  they  then  compute  all  phase  velocities  down  to  whatever 
minimum  period  is  desired.  As  the  second  generation  software  is 
developed,  computations  for  the  group  velocity,  phase  attenuation, 
amplitude  excitation  function,  and  apparent  initial  phase  are  in¬ 
corporated  into  the  procedure  (but  not  yet  into  a  single,  automatic 
routine  combined  with  the  phase  velocity  computations) . 

In  a  series  of  five  later  papers  under  this  contract  (Kausel, 
Schwab  and  Mantovani,  1977;  Mantovani  et  al  1976,  1977;  Mantovani, 
1977  a,b)  this  second  generation  software  was  fully  developed  and 
applied  to  the  generation  of  multimode  theoretical  siesmograms 
containing  as  many  as  21  modes;  each  of  these  was  represented 
over  the  entire  period  range  down  to  1  second. 

The  final  stage  of  our  work  on  the  generation  of  complete 
theoretical  seismograms  for  torsional  waves,  was  the  development 
of  a  third  generation  program.  This  routine  is  fully  automatic, 
and  requires  only  the  structural  specification  as  input.  The  out¬ 
put,  which  is  obtained  in  a  single,  relatively  short  computer  run, 
contains  all  of  the  frequency-domain  information  required 
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to  compute  theoretical  seismograms  for  arbitrary  source  specifi¬ 
cations.  As  desired,  all  body-wave  and  surface-wave  arrivals, 
for  periods  greater  than  ten  seconds,  are  obtained  from  the  90-100 
modes  thus  specified.  Results  from  this  work  were  recently  pre¬ 
pared  and  a  preprint  is  appended  (Liao,  Schwab  and  Mantovani,  1977). 
We  are  now  in  a  position  of  being  able  to  compare  directly,  the 
entire  experimental,  torsional-wave  seismograms  from  the  long- 
period  WWSSN  instruments  with  those  computed  from  theoretically 
specified  sources  and  structures. 

Our  work  on  the  algorithm  and  the  programming  on  the  Rayleigh-, 
or  spheroidal-wave  theoretical  seismograms,  began  with  detailed 
analysis  and  improvements  of  the  basic  direct  method  for  handling 
such  calculations  on  a  sph-rical,  gravitating  earth.  In  its 
original  form,  this  method  was  initially  develoepd  in  the  series 
of  papers  by  Hoskins  (1920) ,  Pekeris  and  Jarosch  (1958) ,  and 
Alterman,  Jarosch  and  Pekeris  (1959);  some  numerical  details  con¬ 
cerning  such  computations  were  given  by  Bolt  and  Dorman  (1961), 
and  later,  by  Takeuchi  and  Saito  (1972).  In  the  appended  manu¬ 
script  (Schwab  et  al ,  1977)  we  describe:  (1)  the  newly  developed 

simplifications  of  the  usual  algorithm,  which  has  made  it  possible  — 
we  believe  for  the  first  time  —  to  develop  a  direct  algorithm  for 
group  velocity  computation,  that  is  independent  of  the  usual 
appeal  to  variational  techniques;  (2)  the  numerical  problems  that 
are  associated  with  this  type  of  computation,  in  quite  detailed 
form,  in  relation  to  what  has  appeared  previously  in  the  literature; 
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(3)  our  optimization  of  the  improved  algorithm,  and  the 
efficiency  of  the  optimization  relative  to  various  similar 
computations  problems. 

Relative  to  the  efficiency  determination,  the  most  per¬ 
tinent  results  affecting  the  work  on  this  contract  are:  (1)  that 
such  computations  —  for  Rayleigh-,  or  spheroidal  waves  on  a 
spherical,  gravitating  earth  —  are  abrut  six  times  more  expensive 
than  the  comparable  torsional-wave  computations,  and  (2)  that 
spheroidal  waves  can  be  handled  on  a  non-gravitating  earth  for  only 
about  twice  the  expense  of  torsional  waves.  From  these  results  we 
conclude  that  the  present  optimization  is  still  too  expensive 
to  use  for  the  computation  of  theoretical  seismograms  for 
spheroidal  waves  (down  to  periods  of  only  ten  seconds) ,  but  that 
the  optimization  technique  will  be  satisfactory  for  this  purpose 
if  a  means  can  be  devised  to  approximate  the  removal  of  gravity 
from  the  formulation.  Such  a  technique  has  already  been  devised 
(Schwab,  1977),  but  the  numerical  tests  have  only  just  begun. 

The  final,  practical  purpose  of  our  work  under  this  contract  — 
application  of  our  results  to  the  discrimination  problem  —  will 
involve  comparison  of  complete  theoretical  and  experimental 
seismograms.  It  is  theorefore  important  that  we  have  as  accurate 
a  means  as  possible  of  obtaining  the  instrumental  constants  from 
the  impulse  response  of  the  experimental  record.  These  constants 
then  permit  us  to  include,  with  a  minimum  of  error,  the  effect 
of  the  instrumental  response  on  the  theoretical  seismogram.  Our 
improved  scheme  for  inversion  of  the  impulse  response  to  obtain 
the  instrument  parameters,  is  described  in  the  appended  preprint 
by  Mitchel  (1977). 
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ABSTRACT 


Algorithmic  and  numerical  analyses  are  carried  out 
for  Ray le igh-wa ve  dispersion  computations  on  a  spherical, 
gravitating  earth.  Our  work  is  based  on  the  direct,  Alterman- 
J ar osch-Peke r is  formulation.  For  practical  purposes,  we  fix 
period  and  determine  the  associated  phase  velocity  (or  polar 
order  number).  Neither  this,  nor  integration  downward  from  the 
free  surface  --  both  "non-standard"  procedures  --  results  in 
unexpected  difficulties.  The  latter  procedure  yields  a  simpli¬ 
fication  of  the  computational  algorithm,  the  clarity  of  which 
al  ws  it  to  be  extended  to  group-velocity  evaulations.  The  AJP, 
d i r ec t - i n t eg ra t ion  formulation  is  optimized  and  compared  with  the 
fastest  --  Knopoff's  method  --  of  the  techniques  based  on  the  flat, 
homogeneous- layer  approximation.  The  optimized  form  of  the  AJP 
method  (spherical)  is  three  times  slower  than  Knopoff's  (flat, 
non-gravitating)  method  when  gravity  is  included  in  the  AJP 
formulation;  and  is  1.36  times  slower  when  gravity  is  not  included. 
Additional  programming  would  reduce  the  former  estimate  to  a  lower 
bound  of  2.42  times  slower,  and  the  latter,  to  a  lower  bound  of  1.30. 
In  size  and  number,  the  treatment  of  integration  "steps"  in  the 
d i r ec t - 1 n t eg ra t i o n  procedure,  is  equivalent  to  the  treatment  of 
"layers"  in  the  homo ge ne ou s -  1 ay e r  approximation;  thus  the  usual 
assumption  that  the  former  method  does  a  better  job  of  treating 
continuous  parameter-depth  distributions,  appears  to  be  invalid. 
Overflow  problems  in  the  AJP  formulation  can  be  controlled  by 


k 


simple  normalization.  Loss-of-precision  problems  appear  to  be 
intrinsic  to  the  AJP  formulation.  At  a  fixed  period,  this  results 
in  the  attainable  accuracy  of  the  phase  velocity  decreasing  as 
mode  number  increases;  and,  for  fixed  accuracy  in  the  phase 
velocity,  as  period  decreases  the  maximum  mode  number  that 
can  be  treated  successfully  decreases. 


1 .  INTRODUCTION 

Dorman,  Ewing  and  Oliver  (1960)  described  Che  use  of  an  elec¬ 
tronic  computer  to  calculate  surface-wave  dispersion  for  multilayered, 
p e r f e c 1 1 y - e 1  a s t i c  half-spaces.  Their  computations  were  based  on 
the  technique  devised  by  Thomson  (1950)  and  Haskell  (1953).  Press, 
Harkrider  and  Seafeldt  (1961)  also  used  the  Thomson-Haskell  tech¬ 
nique,  and  with  a  more  advanced  computer,  greatly  improved  the 
speed  of  computation.  Randall  (1967)  later  applied  Knopoff's 
(1964)  method  to  this  problem  and  reported  a  further  improvement 
in  speed  for  the  Rayleigh-wave  case. 

In  a  later  series  of  papers,  Schwab  (1970)  and  Schwab  and 
Knopoff  (1970,  1971,  1972,  1973)  improved  the  optimization,  for 

computer  application,  of  both  the  Thomson-Haskell  and  Knopoff's 
methods  for  flat,  multilayered  media.  These  papers  also  provide 
complete  details  for  obtaining  full  control  over  the  accuracy  of 
the  computations,  and  for  generalizing  the  algorithms  to  include 
computation  of  attenuation  due  to  the  intrinsic  anelasticity  of 
the  earth. 

For  Love  waves,  the  use  of  spherica.-to-flat  structure  trans¬ 
formations  (Biswas  and  Knopoff,  1970;  Schwab  and  Knopoff,  1971; 

1972;  1973;  Kausel  and  Schwab,  1973)  makes  it  possible  to  carry 

out  all  spherical  dispersion,  attenuation,  and  focal-response 
problems  using  the  optimized  algorithms  for  flat  structures.  Sev¬ 
eral  attempts  have  been  made  to  develope  the  same  type  of  transfor¬ 
mation  for  Rayleigh-wave  computations  (Alterman,  Jarosch  and 
Pekeris,  1961;  Bolt  and  Dorman,  1961;  Biswas,  1972;  Schwab  and 
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and  Knopoff,  1972),  but  these  have  all  yielded  only  empirical  re¬ 
sults  which  lack  general  applicability.  Thus,  at  the  present  time 
at  least,  it  appears  that  one  cannot  apply  transformation  theory  to 
Rayleigh -wave  dispersion  computations  on  any  arbitrary,  spherical, 
gravitating  earth.  Bhattacharya's(1976)  recent  results  --  although 
we  will  not  pursue  this  approach  in  the  present  paper  --  suggest  the 
feasibility  of  an  interesting  new  procedure  for  treating  spherical, 
gravitating  structures:  Gravitation  alone  might  be  handled  by  trans¬ 
formation  techniques,  while  Bhat tacharya ' s  approach  could  be  used  to 
optimize  the  treatment  of  sphericity. 

The  primary  purposes  of  the  present  paper  are:  (1)  to  report 

on  our  study  of  the  optimization  of  the  direct  computations  (see 
Wiggins  (1976)  for  a  discussion  of  computations  based  on  the 
variat ional  technique),  (2)  to  report  the  results  of  our  study  con¬ 
cerning  accuracy  considerations,  and  (3)  to  determine  the  efficiency 
of  these  direct  computations  relative  to  the  analogous  computations 
for  non-gravitating  structures.  Also,  a  new  computational  technique 
is  developed,  for  the  calculation  of  group  velocities,  which  does 
not  depend  on  the  numerical  evaluation  of  "energy  integrals."  Our 
second  purpose  is  to  present  --  we  believe  for  the  first  time  --  an 
explicit,  quantitative  comparison  of  the  relative  efficiencies  of 
the  two  basic  techniques  for  performing  surface-wave  dispersion 
computations:  that  in  which  an  exact  structural  specification  is 

employed  with  approximate  mathematical  methods,  and  that  in  which 
exact  analytical  techniques  are  applied  to  an  approximate  model  of 
the  structure,  i.e.,  where  the  structure  is  replaced  by  a  sequence 
of  homogeneous  layers. 


ALTERMAN-JAROSCH-PEKERIS  (AJP)  FORMULATION 
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The  basic  formulation  for  our  problem  (Pekeris 
and  Jarosch,  1958)  involves  the  solution  of  three  s econd -o rder , 
ordinary  differential  equations  constrained  by  a  set  of  boundary 
conditions.  For  purposes  of  numerical  solution  it  is  advisable 
to  reduce  this  system  to  six,  linear,  first-order  differential 
equations,  as  was  done  by  Alterman,  Jarosch  and  Pekeris  (1959). 

Bolt  and  Dorman  (1961)  applied  this  formulation,  to  the  evaluation 
of  Ray 1 e i gh -wave  dispersion,  and  reported  on  those  numerical 
details  which  it  was  economically  feasible  to  investigate  with 
second-generation  computing  equipment.  Detailed  algorithmic  testing 
of  accuracy,  precision,  and  efficiency  characteristics  really 
requires  the  present,  third-generation  machinery,  which  we  have 
employed  in  the  current  study;  the  work  we  report  here  can  be 
considered  as  the  logical  extension  of  the  above  series  of  papers. 


To  sketch  the  AJP  formulation,  if  we  let 
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with  y;  and  y3  related  to  the  components  of  displace¬ 
ment  u  (r  ,  6, *) ,u  (r  ,  6  ,  $) ,  and  u .  (v,8,$)  by 

r  «  <p 

i  \  iu  t 

u  =y i (r)x  ,{0)e  e 
r  t 


,  \  d  m  /  o  %  imp  i-w  t 
%  =  y3(r)d0xf  (6)e  e 


(2.02) 


im  ,  ,  ...  im$  iut 

V  77^  y3(r)xa(9)e  e 


For  propagating  surface  waves  diverging  from  the  epicenter. 


Xm=  i( pm+i— Qm) , 

K  i.  2  ft  n  ' 


(2.03) 


for  waves  converging  toward  the  epicenter, 

m  1  ,  m  .  2„m. 

X  “tCP.-x-Q.) • 

11  i  it  l 


(2.04) 


For  a  treatment  of  the  situations  which  require 
the  use  of  (2.03),  (2.04),  or  their  sum,  see  Schwab  and  Kausel 

(1976).  In  this  same  reference,  the  justification  is  given  for  our 
major  departure  from  previously  reported  computations  of  Rayleigh 
wave  dispersion  on  a  spherei  Strictly  speaking,  Rayleigh  waves 
only  exist  on  a  sphere  at  the  discrete  set  of  frequencies  corres¬ 
ponding  to  integral  values  of  the  polar  order  number  _2  .  However, 
fixing  and  evaluating  the  corresponding  angular  frequency  cj  does 
not  yield  the  dispersion  data  at  equal  frequency  intervals,  which 
we  desire  to  use  in  the  usual  numerical  technique  for  obtaining 
time  series  by  inverse  Fourier  transformation.  Schwab  and  Kausel 
(1976)  have  shown  that,  for  most  practical  applications  of  propaga¬ 
ting  surface  waves,  non-integral  l_  at  equally-spaced  frequencies 


can  be  used  without  introducing  significant  errors;  therefore, 
we  adopt  the  procedure  of  fixing  m  and  computing  £,  or 

c  =  aw/(l+l/2) ,  (2.05) 

where  a  is  the  radius  of  the  earth.  The  relation  between  £ 
and  the  true  spherical  phase  velocity  is  also  treated  by  Schwab 
and  Kausel  (1976).  In  equation  (2.01),y2  and  y4  are,  respectively, 

the  radial  dependences  of  the  r_r,  and  the  rQ_  and  r_£  components  of 
stress;  y  and  y6  arise  from  the  presence  of  tne  gravitational 

field. 

Since,  in  any  numerical  integration  procedure, 
it  is  important  to  initiate  the  integration  with  accurate  values, 
we  have  chosen  to  proceed  from  the  free  surface  downward.  This 
allows  us  to  specify  the  initial  vector  exactly.  The  integration 
is  then  carried  down  to  a  depth  sufficient  to  make  it  immaterial — 

to  n  significant  figures  in  £  or  £ - just  how  we  terminate  the 

integration:  with  an  approximation  of  a  free  surface  or  rigid 

surface,  for  example.  The  fact  that  such  a  termination  process  is 
valid  has  been  checked  by  extensive  numerical  tests  in  the  course 
of  this  work.  These  tests  follow  the  lines  of  the  1 a y e r - r ed uc t  i  o n 
experiments  described  by  Schwab  and  Knopoff  (1970;  1972),  and  will  be 

described  in  some  detail  below.  In  Section  7  we  discuss  the  termina¬ 
tion  of  the  structure  at  depth  by  either  a  solid  or  liquid,  homogen¬ 
eous,  gravitating  sphere. 


Here,  we  should  point  out  that  the  warning 
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given  by  Takeuchi  and  Saito  (page  24  1  ,  1  9 7 2 )  '  aga i ns t  proceeding 

downward  from  the  free  surface  when  integrating  the  system  of 
differential  equations,  or  when  forming  the  layer-matrix  product 
if  applying  the  Thomson-Haskell  technique  or  Knopoff's  method, 
does  not  appear  to  be  justified  by  our  experience.  In  the  work 
upon  which  we  report  herein,  downward  integration  did  not  give 
rise  to  any  unexpected  difficulties;  in  previous,  extensive  work 
witli  matrix  methods  applied  to  Rayleigh-wave  dispersion  computations 
(Schwab,  1970;  Schwab  and  Knopoff,  1970;  1972),  the  formation  of 
matrix  products  upward  toward  the  free  surface  (Thomson-Haskell 
formulation)  was  not  found  to  have  any  advantage  over  formation  of 
the  product  in  the  downward  direction  (Knopoff's  formulation). 


y  vanish , 

_ U 


and 


Continental  structure.  In  this  case,  y  and 
■- y  (  £+1 ) / a  at r  *  a.  Thus  we  can  write  the  starting 


vector  as 


I 


o  r 

V  (a)*sy1  (a)Xi  (a)+y3(a)X2(a)+y5(a)X3(a)  ,  (2.07) 

arid  for  £■'  a 

Vs (r ) «y i (a) Xi ( r ) +y 3 (a) X2 (r)+y5 (a)X3 ( r)  .  (2.08) 

The  three  quantities  which  are  unknown  --  yi(a),  yj  (a)  , 
and  y^fa)  --  can  be  carried  implicitly  in  the  computations, 

while  we  integrate  the  vectors  whose  starting  values  are  known 
exactly:  X;,  X2 ,  and  X3.  That  is,  we 

integrate  to  obtain  Xi  at  depth;  this  is  repeated,  in  turn, with 
X2  and  X3.  Thus  we  actually  use  e qua t io n ( 2 . 0  1 )  in  the  form 
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X  =AX  to  integrate  from  the  surface  r  «  a  to  the  depth  at 
i  i  —  — 

which  the  boundary  conditions  are  to  be  applied:  r  =  r^  ,  where 

we  can  again  express  Y<.  in  terms  of  the  undetermined  coefficients 

by  using  equation  (2.08). 

If  we  define  a  rigid  boundary  at  depth  by 
yi  (r0)=y3(ro)=y!>(ro)=0,  (2.09) 

we  then  obtain  three  linear,  homogeneous  equations  in  three 

unknowns - the  undetermined  coefficients - and  the  determinant  of 

the  coefficient  matrix  must  vanish  if  we  are  to  have  a  non-trivial 
solution.  Thus  the  dispersion  function,  F  ,  takes  the  form 

[X 1  (  r  o  )  ]  i  [x2U0)]i  [x3Uo)Jl 

Fa(c,u»)  -  [x  i  ( r  o )  ]  3  [X2(r0)i3  [X3(r0)]j  ,  (2.10) 

[X  i  ( r  o  )  ]  5  [X  2  (  r  0  >  1  5  [X  3  (  r  o  )  ]  3 

zeros  of  which  define  valid  (cfU))  dispersion  pairs.  For  the  two 
approximations  to  free  boundaries  at  depth,  we  have  used  the 
definitions 

Y2 (Co)  “  Vu  (to)  “  0  (2.11) 

and  either 


y 6  ( r o )  y  s  (r  o)  +  D  /r  o 


(2.12) 
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o  r 


y 6<r, )=-ys(r0) U+l) /a  .  (2.13) 


which  yield,  respectively,  dispersion  functions 
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Oceanic  Strut  mrc .  In  this  case,  the  analog  of  equation 
(2.01)  is,  for  the  homogeneous  oceanic  (liquid)  layer. 
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At  £  ■  a,  y2  vanishes  and  y 6 “_y 5 (£+1 ) /a ,  and  we  can  write  the 


starting  vector  as 
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YL(a)  =  y ] (a)  Zj(a)  +  y5(a)  Z2Ca),  (2.18) 

and  for  r<a 

Y  ^  (r )  =  y  i  (a)  Z^r)  +  y5(a)  Z2(r).  (2.19) 


Again,  we  carry  the  unknown  quantities 


y,  (a)  and  y5  (a)-- 


implicity,  and  integrate  the  vectors  whose  starting  values  we 
know  exactly:  and  Z2,  using  equation  (2.16)  in  the  form  Z^  =  B  Z^. 


On  the  oceanic  side  of  the  liquid-solid  boundary  at  the  bottom 


YL(r, )=y  j (a) 


_r  *  r  i  ,  we 

then  have 

[z  i  ( r  i  >  j  i 

[z  2  (  r  1  )  ]  1 

[Z  1  (  r  1  )  j  2 

+  y 5 (a) 

lz2  (  r  1  )  1  2 

lzl  Cr  j  )  j  5 

Lz  2  (  r  1  >  1  S 

[Zl(ri)]6 

Cz  2  (  r  ]  )  ]  6 

■-  - 

_  _ 

(2.20) 


At  this  boundary,  ylf  y2,  ys,  anu  y^  are  continuous,  yi 4  must 


vanish,  and  y3  is  undetermined.  Thus  on  the  solid  side  of  this 
interface,  we  have 
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3.  ALGORITHM  FOR  GROUP  VELOCITY  DETERMINATION 


When  treating  Rayleigh  waves  on  a  spherical,  gravitating 
earth,  the  variational  technique  is  usually  employed  to  compute 
group  velocities  (Takeuchi  and  Saito,  1972,  Section  III).  This 
involves  the  evaluation,  by  numerical  approximation,  of  integrals 
having  the  form 

■a 

f ( r )  yi(r)y.(r)  dr  .  (3.01) 

I  J 

Jo 

Since  the  functions  yk(r)  become  highly  oscillatory  for  large 
(radial)  mode  numbers,  this  numerical  evaluation  can  become 
inaccurate  (Knopoff  £t^  a_l  .  ,  1974,  Appendix).  Further,  I  y^ ( r  )  i 

become  spuriously  large  at  depths  much  below  those  at  which  there 
is  significant  energy  in  the  mode,  at  the  period  being  treated. 
Although,  somewhat  surprisingly,  these  spurious  magnitudes  do  not 
seem  to  affect  the  location  of  a  root  of  the  dispersion  function, 
they  can  cause  large  errors  in  the  evaluation  of  integrals  such  as 
(3.01).  Thus  one  must  specify  r2  ,  the  value  of  r  below  which 
the  spurious  magnitudes  occur,  prior  to  seeking  the  group  velocity, 

and  then  evaluate 

f(r)yi(r)  y^r)  dr  (3.02) 

in  place  of  (3.01).  However,  without  prior  knowledge  of  the  group 
velocity,  can  be  quite  difficult  to  determine  in  period  ranges 
such  as  those  in  which  the  energy  shifts  back  and  forth  repeatedly 
between  the  crustal  wave  guide  and  the  low-velocity  channel  in  the 
upper  mantle  (Panza,  Schwab  and  Knopoff,  1972;  Frantsuzova,  Levshin 
and  Shkad inskaya ,  1972;  Schwab  and  Knopoff,  1971;  1972).  Since 

the  integrals  must  eventually  be  evaluated  to  obtain  excitation 
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and 


•  £>  p 

F  -  <TF> 


(3.06) 


is  given  by  the  same  type  of  expression  as  (3.05),  with  dots 
replacing  the  primes.  The  elements  |X^(r0) in  (3.05)  and 
(3.06)  are  obtained  exactly  as  described  in  Section  2.  The 
evaluation  of  i X !  ( r o )  | ^  and  |X^(r0)j  requires  a  simple  extension  of 

the  algorithm. 


C  o  nt  inrnt.il  Structure  .  In  this  case  we  start  with  the 
s i x t  h - o  r  d  e r  system 


and  form 


Y=  A  Y  , 


Y '  =  A ' Y+AY  ’ 

*7  *  • 

Y  =  AY  +  AY. 

Here  again,  we  sn  these  equations  of  motion  in  terms  of  the 

X  ,  X’  ,  X  .  t!  i  we  know  exactly  at  the  surface: 
i  l  i 


(3.07) 

(3.08) 

(3.09) 
vectors  , 


A’  X .  +  AX’  . 

l  i 


(3.10) 


I.  =  A  X .  +  AX . 
1  11 


(3.11) 


Since  X.  can  be  determined  i  nd  e  p  e  i  1  e  n  1 1  y ,  we  car.  treat  A'X^  and 


x  as  known  vectors  at  each  depth,  and  we  have 
i 


L 


X'  =  AX'  .  +  C  . 
i  11 


X  .  =  AX .  +  D . 

1  li 


where 


C.  ( r )  =A '  (r)X.(r)-2p(r), 


[Xi<r)Ji 

0 

[X.(r)j3 


D  (r)=AU)X.  (r  )  =-  §4- 
i  i  c 4  r 


(2e  +  D 


[>./(A+2i)]  [X.(r)j3 

[og-2u(3X+2u)/(A+2y)r]  [> 


[4u(X+u)/ (A+2u)r]  [X.(r)j 

0 

[-4nG;-j  [x  (  r  )  ]  3+ [l  /  r]  [x 
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(3.12) 

(3.13) 

(3  .14) 

1. (r)] 3+[Xj (r) 

3 

(r)]  5 


(3.15) 


(3.19) 


which  can  be  written  in  terms  of  known  vectors,  B’Z.  and  BZ  .  , 

1  1 


at  each  depth : 


Z'.  *  BZ\  +  E.  (3.20) 

l  li 


Z .  »  BZ.  +  G.  ,  (3.21) 

l  l  l 


whe  re 


E.(r)=B'(r)Z.(r)=[2i(i+l)/r2u3j 


-H . (r) 
1 


Ii(r) 


[4nGp(r>]  H.(r) 


(3.22) 


H.  (r) 


G.(r)=B(r)Z^(r)  =  [-a(2f.+  l)/c2r2uj] 


[g  (r)p  (r)[Hi  (r) 
0 


JjL(r) 


(3.23) 


and 


Hi ( r) -g (r)  [Zi(r)j j  -  [l  /  o  (  r  )]  [zi(r)]2 


[Z±  <^)]  5 


(3.24) 
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vectors,  i.e.  approximately  twice  the  amount  of  time  it  takes  to 
compute  a  single  value  of  the  dispersion  function  while  iterating 

for  £ .  The  number  of  iterations  required  in  the  computation  of  £ 
is  variable,  but  for  the  present  purpose,  ten  iterations  may  be 

taken  as  a  representative  value.  Thus  30  integrations  would  be 

required  to  obtain  c,  and  only  6  to  obtain  u.  To  the  accuracy  of 
this  estimate  then,  £  can  be  obtained  in  only  20  percent  of  the  time 

the  evaluation  of  £  requires. 
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4 .  NUMERICAL  TECHNIQUE  FOR  INTEGRATING  THE  SYSTEM  OF  DIFFERENTIAL 
EQUAL  IONS 

To  optimize  the  computations,  it  is  useful  to  constrain 
the  structural  specification  somewhat.  For  this  purpose: 

1.  The  liquid,  oceanic  layer  is  limited  to  a  single, 
homogeneous  layer.  A  special  fourth-order  Runge-Kutta  technique 
(see  below)  is  used  for  the  first  three  steps  of  the  integration-- 
step  size  of  about  I  km--  and  a  fourth-order  pred i cat or-cor rec t o r  method  (see 
below)  is  employed,  if  necessary,  to  continue  the  integration  to 
the  bottom  of  the  oceanic  layer  with  the  same  step  size. 

2.  The  sedimentary  layers  are  limited  to  a  sequence  of 
homogeneous  layers,  each  of  which  does  not  exceed  1  km  in  thickness. 
These  layers  are  treated  with  one  fourth-order  Runge-Kutta  step. 

3 .  The  sub  sedimentary  crustal  layers  must  also 

be  homogeneous,  and  each  is  treated  as  tiie  oceanic  layer;  the  step 
size  fixed  at  about  1  km. 

4.  The  sub-Moho  mantle  requires  continuous  velocity- 
depth  and  dens ity-depth distributions  ,  although  discontinuities 


l 


can  be  approximated  as  closely  as  desired  by  specifying  large 
gradients.  As  with  the  oceanic  layer, three  Runge-Kutta  steps 
are  followed  by  a  fourth-order  predictor-corrector  method.  The 
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initial  step  size  is  1.5625  km,  with  which  we  execute  three  Runge- 
Kutta,  and  seven  predictor-corrector  steps.  The  step  size  is  then 
doubled  andfive  p r ed ic t o r -co r r ec t o r  steps  are  perft  'med ;  this 
procedure  is  repeated  until  the  step  size  reaches  12.5  km,  and  the 
predictor-corrector  method  is  then  applied  with  this  fixed  step 
size.  The  results  of  our  numerical  testing  have  shown  that  this 
dependence  of  step  size  on  depth  is  sufficient  to  yield  4-signif_i 
cant-figure  accuracy  in  the  computed  values  of  c.  Concerning  this 
point,  one  should  review  the  treatment  given  by  Schwab  and  Knopoff 
(1972),  in  which  piecewise-cont inuous  velocity-  and  density-depth 
distributions  are  treated  with  the  homo geneo u s -1  ay e r  approximation. 
Comparison  will  show  that  the  thickness  of  the  layers  as  a  function 
of  depth,  in  that  approximation,  is  roughly  the  same  as  the  integration 
step  size  in  the  present  analysis.  That  is,  to  the  degree  of  accuracy 
which  two  such  dissimilar  methods  can  be  compared,  if  the  step  sizes  above 
are  used  as  layer  thicknesses  in  the  homogeneous  layer  approximation,  that 
technique  will  yield  4  significant  figures  in  the  computed  values  of  c_  . 
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Runge-Kutta  technique.  To  start  the  predictor-corrector 
method  we  have  used  a  Runge-Kutta  technique  designed  specifically 
for  this  purpose.  Here  one  is  only  interested  in  being  able  to 
minimize  the  hounds  on  the  truncation  error.  Ralston  (1962)  has 
treated  this  problem,  and  gives  the  following  algorithm  for 
o  h  t  .lining  the  first  four  points  for  starting  a  predictor- 
corrector  method. 


In  terms  of  a  single,  first-order  differential  equation; 


y  =  £  (r,y)  ,  y(rQ)  =  yQ  , 

at  r,,  r2,  ....  the  Runge-Kutta  method  is  given  by 

yn+ryn  Vf 


Here,  y  =  y(r  ),the  w.  are  constants,  and 
n _ n  i 


i  - 1 


K.=h  f(r  +  a . h  , y  +  .  1 ,  6  .  .  K  ) 
i  n  n  in  ti  j  =  l  ij  j 


with  h  =r  -r  ,  and 
n  n  + 1  n 


a  i  =0 
a 2  =  2/5 

u  3  =  (14-3/T)/16 

a  4  =  1 


(4.1) 

(4.2) 


(4  .3) 


(4.4) 


% 
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32  1 

=  a  2 

S31 

=  a  3  - 

6  32 

3  32 

=  [a  3 

(a  3-a2)]/  [2a2 

( 1 -2a  2  )  ] 

341 

=  a4- 

3  4  2  —  3  4  3 

-  (1- 

aa)ra->+a3  -1- 

(2a, -1)21 

&  4  2 

2a2 

(a  3-a2 )  [6a2a  3 

—  4  ( a  2  +ci  3  )  +  3  j 

64  3 

(1- 

2<»g)  (l-a2)  (1- 

a3) 

013  (a3-a2)  [6a2a3 

-4(a2+a3)+3] 

wl  't  +  [l -2  ( a 2+a  3)^ 2  a2a3 

w2  =  (2  a  3  -1) / [l2a2(a3-a2) (l-a2)] 

(4.6) 

3  =  (1-2  a2)/[l2a3  (  a  3~  o2)  (  1 -a  3  )  ] 

*  _i_  +[2(a2+a3)_3j/[l2(l-a2)  (l-cjj)] 

Predictor-corrector  method,  jhe  fourth-order  method  we 
have  used  (Hamming,  1959)  is  fully  described,  along  wich  the  details 
concerning  doubling  of  step  size,  by  Ralston  (1960).  Highly  practi¬ 
cal  details  concerning  the  combination  of  the  Runge-Kutta  and  predic¬ 
tor-corrector  routines,  which  we  have  employed,  will  be  found  in 
Anon. (1970).  One  should  be  warned,  however,  that  the  use  of  the 
subroutines  therein  described  is  highly  inadvisable  for  our  present 
purposes.  The  use  of  these  general  purpose  subroutines  can  increase 
computation  expense  by  a  factor  of  10  to  100  over  that  of  the  optimized 
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OPTIMIZATION  OF  THE  ALT EKMAN- J AROS CH- PEKEKI S  FORMULATION 
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The  key  to  optimizing  the  integration  is  to  apply  our 

knowledge  about  this  specific  problem  to  specify  all  the  depths,  r^> 

at  which  a..(r)  are  to  be  evaluated.  The  evaluation  of  these 

U 

elements  can  then  be  removed  from  the  innermost,  integration  loops 

of  the  program.  The  details  concerning  these  depths  are  contained 

in  the  preceding  section.  In  Figure  1,  the  optimized  scheme  for 

the  evaluation  of  a..(r  )  —  for  the  solid  sedimentary  layers,  subsedi_ 
1 J  k 

meniarv  crustal  layers,  and  mantle — is  indicated  in  outline  form. 

This  figure  shows  that  most  of  the  procedure  for  evaluating  a..(r.) 

— l.i _ _ 

can  even  be  removed  from  within  the  wand  £  loops:  within  the  £ 

loop,  each  new  £  value  requires  only  bli+l  assignments,  6N  +  1  multiply 

cations,  and  N+l  additions;  within  the  w  loop,  each  new  w  value 

requires  only  3N  +  1  assignments,  N+l  additions,  and  2N^  subtraction^, 

where  N  is  the  number  of  depths  at  which  a  .  .  (  r  )  must  be  evaluated. 

1 J  k 

All  other  portions  of  the  element  determinations  are  performed 

external  to  these  loops.  In  Figure  2  the  same  information  is  given 

for  the  elements,  b. ,(r)  of  the  matrix  describing  the  integration 

_JJ _ 

through  the  liquid,  oceanic  layer.  Again,  all  elements  can  be 
evaluated  external  to  the  integration  routine. 


In  the  integration  procedures  themselves,  it  is  very 
important  to  form  matrix  products,  such  as  those  in  (2.01)  and 
(2.16),  in  an  explicit  manner.  This  permits  full  use  to  be  made 
of  the  many  zero  elements,  and  those  that  are  independent  of  £, 
or  are  equal  to  another  matrix  element.  For  example,  the  basic 
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AJP  matrix  multiplication  for  solid  layers  is  illustrated  in 
Figure  3. 


In  reports  on  our  earlier  work  with  Love-wave  dispersion 
and  dispersion-attenuation  computations,  for  both  flat  and 
spherical  structures,  it  was  possible  to  give  the  short,  key, 
FORTRAN  program  segments  (Figure  2,  Schwab  and  Knopoff,  1972; 
Figures  4  and  5,  Schwab  and  Knopoff,  1973).  These  together 
with  descriptions  of  the  root-bracketing  and  root-refining  proce¬ 
dures,  completely  specify  the  optimization  when  the  mul t i- , homoge¬ 
neous-layer  approximation  is  employed.  When  this  approximation  is 
used  with  Rayleigh  waves  on  flat  structures,  the  optimization  can 
be  specified  in  the  same  manner  (Figures  11,  12  and  13,  Schwab  and 
Knopoff,  1972).  When  employing  the  method  of  direct  integration 
of  the  equations  of  motion,  it  is  not  possible  to  exhibit  the 
complete  program  optimization  in  this  compact,  simple  manner  for 
Rayleigh-wave  dispersion  on  a  spherical,  radially  heterogeneous, 
gravitating  earth.  However,  it  i^s  possible  to  present  the  most 
important  part  of  the  algorithm  as  a  relatively  compact  program 
segment.  This  is  given  in  Figure  4a, which  illustrates  the  predic¬ 
tor-corrector  method  we  have  applied  to  the  integration  from  below 
the  Moho  to  the  selected  value  of  rQ;  the  automatic  doubling  of 
integration  step  size  is  included  in  the  segment.  Most  of  the 
computation  time  is  spent  in  this  program  segment,  which  is  entered 
with : 

C0EFF1  -  (4/3)  H 
C0EFF2  =  311 

C0EFF6  =  (121/36)  H 
H  -  -25/16; 

the  indices  for  the  successive  integration  step-size  regions  are 
given  in  Table  1;  for  the  details  concerning  B(I,J),  see  the 
description  of  subroutine  DHPCL  (Anon.,  1970).  In  this  type  of 

Programming  there  are  important,  machine-dependent  considerations: 


the  manner  in  which  the  6x6  matrix  elements  are  stored  in 
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memory,  and  the  manner  in  which  indices  are 

handled  in  the  program  segment  given  in  Figure  4a.  In  fact, 
the  indices  ITP1,  ITP3,  ITP8,  ITP9,  ITP10  which  are  used  for 
compactness  in  Figure  4,  actually  slow  computation  speed  on 
the  IBM  360/91;  these  subscripts  are  best  used  in  explicit 
form  IT+1,  IT +3,  I T+8 ,  IT+9,  1T+10.  The  key  program  segment  is 

given  in  the  form  shown  in  Figure  Aa  for  two  reasons:  (1) 

to  illustrate  the  logic  as  clearly  and  simply  as  possible, 
and  (2)  to  provide  an  illustrative  example  of  the  importance 
of  handling  subscripting  and  storage  in  the  manner  most 
appropriate  for  a  given  machine.  The  time  required  to 
execute  DO-loop  170  once,  specifies  the  necessary  time  to 
execute  one  integration  step  for  each  of  the  three  vectors 
X.  ;  thus,  to  perform  one  integration  step  in  forming  the 
dispersion  function,  DO-loop  170  must  be  executed  three 
times.  The  time  for  one  integration  step  in  forming  F_  is 
termed  the  characteristic  time  2.»  which  we  use  to  illustrate 
t lie  importance  of  correct  subscripting  and  storage.  The 
characteristic  time  for  the  segment  in  Figure  Aa  is 
A89  nsec /s t ep / i t e r a t ion .  By  simply  reversing  the  order  of 
the  subscripts  of  B(I,J),  this  time  is  improved  by 
92  lisec/step/iteration;  if  ITP1,  etc.  are  used  explicitly 
as  IT+1,  etc.,  is  decreased  still  further  by  an  amount  of 
67  Msec/step/iteration;  and  if  a^ .  are  stored  more  efficiently, 
still  another  AA  u s ec / s t e p / i t er a t ion  can  be  saved,  bringing 
T  down  to  286  p s e c / s t ep / i t e r a t i on .  DO-loop  170,  in  a  form 
incorporating  the  above  improvements,  is  shown  in  Figure  Ab. 
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It  should  be  understood  that  286  psec/step/iteration  is  a 
lower  bound ;  the  corresponding  effective  cha rac t e r i s t i c  time 
given  in  (5.02)  --  must  reflect  time  spent  in  other  parts  of 

the  program  than  DO-loop  170,  e.g.,  time  spent  in  starting  the 
predictor-corrector  scheme  with  the  Runge-Kutta  procedure. 

Our  computations  have  been  performed  on  IBM  360  and  370 
installations:  a  360/91  in  Los  Angeles,  a  360/65  in  B„ri,  and 

a  370/145  in  Santiago  and  in  Cosenza.  The  first  installation 
was  also  used  in  the  final  optimization  of  the  flat-structure 
Rayleigh-wave  work  (Schwab  and  Knopoff,  1972),  which  allows  us 
to  make  an  accurate  evaluation  of  the  iolative  characteristic 
times  (Schwab  and  Knopoff,  1972)  for  computations  with  flat, 
non-gravitating  structures  and  with  spherical,  gravitating 
models.  In  the  former  case,  this  time  is 

T  =110  usec/layer /iteration  (5.01) 

FLAT  RAYLEIGH 

(which  corresponds  to  Knopoff's  method  applied  to  a  sequence 
of  homogeneous  layers),  and  in  the  latter  case. 


SPHERICAL  T 


RAYLEIGH 


=  336 


psec/step/iteration 


(5.02) 


(which  corresponds  to  the  optimization  of  the  AJP  formulation 
herein  described).  Again,  these  characteristic  times  were 
measured  on  the  IBM  360/91  at  UCLA.  As  we  have  noted  above. 
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an  integration  "step"  can  be  considered  nearly  equivalent  to 
a  "layer"  in  computations  based  on  the  homogeneous- lay er 
approximation.  Also,  to  the  accuracy  possible  in  this  type  of 
comparison,  the  "iterations"  required  in  the  two  cases  (see 
Schwab  and  Knopoff  (1972)  for  details  concerning  iteration 
procedures)  can  be  considered  equivalent.  Thus,  the  relative 
efficiences  of  the  two  types  of  Rayleigh-wave  dispersion 
computation  can  be  evaluated  by  simple  comparison  of  their 
characteristic  times,  and  we  find  that  the  inclusion  of 
sphericity  and  gravitation  triples  actual  computation  time. 

Thus,  the  time  required  to  integrate  each  of  the  three  vectors 

over  depth  in  the  spherical,  gravitating  case,  is  the  same  as 

that  required  to  carry  out  the  analogous  operation  --  the 

formation  of  the  matrix  product  --  for  the  flat,  non-gravitating  case. 

To  obtain  a  valid  comparison  of  the  d i re c t- in t egr a t ion 
method,  with  the  homog e ne ou s- lay e r  technique,  clearly  we  should 
not  include  gravity  in  the  former  method.  The  removal  of 
gravity  reduces  the  vectors  X^,  that  must  be  integrated,  from 
three  to  two,  and  the  number  of  elementary  operations  (multi¬ 
plications  and  additions)  in  (2.01)  from  34  to  23.  Thus,  the 
ratio  of  computation  times  for  non- gr a vi t a t ing  and  gravitating 
spheres,  when  treated  with  the  d i r ec t- in  teg r a t ion  method,  is 
approximately  (2x23)  /  (3x34);  or,  for  the  non-gravitating  case. 
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151  u se c / s t ep / i t e ra t i on . 


(5.03) 


Thus  the  direct  integration  met.iod,  for  a  non-gravitating 
sphere,  is  only  36  percent  slower  than  the  optimized  com¬ 
putations  for  a  flat,  non-g ra vi t a t i ng  structure,  where  the 
latter  is  treated  with  the  homogeneous-layer  approximation. 

From  this  result  we  are  led  to  conclude  that  attempts 
to  devise  Rayleigh-wave  transformation  techniques  --  which 
have  hitherto  been  concentrated  on  curvature  corrections  to 
permit  spherical,  gravitating  structures  to  be  treated  with 
algorithms  for  flat,  non-gravitating  models  --  might  also 
be  directed  toward  gravity  corrections  that  would  allow  one 
to  use  programs  for  spherical,  non-gravitating  structures  to 
handle  the  spherical,  gravitating  case. 

A  final  improvement  in  the  lower  bounds  of  the  A J  P 
characteristic  times  is  possible.  The  limiting  bound  for 
the  gravitating  case  was  obtained  by  including  the  integration 
of  all  three  X^,  simultaneously,  within  DO-loop  170;  the  re¬ 
sult  was  266  p s e c / s t ep / i t e ra t ion .  If  the  two  vectors  of  the 

non-gravitating  case  are  handled  simultaneously  within  the 
loop  (see  Figure  4c),  the  lower  bound  of  Jt  becomes 
143  u sec / s t ep / i t e r a t ion . 

It  is  obviously  of  interest,  to  those  involved  in  surface- 
wave  computations,  to  have  some  idea  of  the  relative  speeds  of 
such  computations  for  the  various  computers  currently  in  use  at 
the  larger  installations.  In  Table  2  we  present  the  results  of 
a  first  effort  to  summarize  this  information  for  the  computers 
that  are  available  to  us  for  testing.  The  test  routine  we 
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have  employed  is  that  given  by  Schwab  and  Knopoff  (Figure  2, 

1972).  A  100-layer  structure  was  used;  the  program  segment 

was  enclosed  within  a  DO-loop  that  was  executed  either  100  or 

1000  times;  case  1  involved  tests  with  c  <  0  ,  and  case  2,  tests 

—  m 

with  c>0  .  Results  from  a  more  extensive  set  of  relative  computer- 
—  m 

time  tests,  for  both  Love-  and  Rayleigh-wave  dispersion  com¬ 
putations  are  given  by  Porter  e_t  al_.  (  19  7  7). 


EXISTENCE  OK  SOLUTIONS  AS  A  FUNCTION  OF  NUMERICAL  AND 


6  . 

ALGORITHMIC  PROCEDURES 

On  IBM  360  equipment,  large-scale  numerical  work  is 
routinely  carried  out  in  double  precision:  about  16  decimal 

digits.  Except  where  indicated  otherwise,  this  precision  was 
used  to  investigate  the  existence  criteria  for  solutions  from 
our  optimization  of  the  basic  AJP  formulation. 

Our  testing  procedure  followed  the  lines  of  the 
1  a y c r - r e d u c t  i  on  experiments  described  by  Schwab  and  Knopoff 
(1972):  At  each  of  a  set  of  periods,  £  is  computed  for  a  com¬ 

plete  range  of  terminating  values,  r0,  for  the  integration.  At 
each,  fixed  period,  by  comparing  the  values  of  c  as  a  function 
of  r_£,  tlie  range  of  r_o_,  over  which  £  is  stable  to  4  significant 
figures  is  immediately  evident.  In  terms  of  _r_^  and  period,  our 
results  for  an  oceanic,  and  a  continental  shield  structure 
(Figure  5)  are  given  in  Figure  6;  the  fundamental  and  first 
seven  higher  modes  are  treated  in  each  case. 

The  results  are  similar  to  those  previously  obtained  for 
Rayleigh  waves  when  the  homogeneous-layer  approximation  is 
employed  (Figures  14,  15,  17,  Schwab  and  Knopoff,  1972):  At 

each  period,  a  certain  minimum  amount  of  structure  (maximum  r  p ) 
must  be  retained  to  ensure  4  significant  figures  in  £.  This 
maximum  value  of  _r_^  is  a  physical  limitation.  For  the  mode  and 
period  of  interest,  there  is  significant  energy  content  down 
to  a  depth  of  a- r o  ,  and  the  structure  above  this  point  thus 
affects  all  4  figures  of  £. 
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In  the  initial  program  testing,  it  is  useful  if  one  can 
integrate  to  any  depth,  irregardless  of  whether  this  results  in  loss 
of  accuracy  in  the  computed  phase  velocity,  or  in  the  expenditure  of 
more  computer  time  than  if  the  integration  were  terminated  at  the 
minimum  acceptable  depth  for  the  desired  accuracy  in  the  phase  velo¬ 
city.  The  use  of  a  very  large  value  of  r^,  or  more  precisely,  a 
large  number  of  wavelengths  of  structure,  will  result  in  overflow; 
thus  a  simple,  temporary  solution  to  this  problem  is  useful.  Such  a 
solution  is  the  simple  extension  of  the  normalization  technique 
described  by  Schwab  and  Knopoff  (1970;  1972).  The  application  of 
normalization  to  the  direct-integration  procedure  is  quite  simple; 
it  is  no  c  necessary  to  begin  normalization  until  the  application  of  the 

predictor-corrector  method  has  begun  in  the  mantle.  To  normalize, 
all  one  need  do  is  determine  the  maximum  of  the  absolute  values  of 
y  ^  at  the  end  of  each  integration  step;  one  then  divides  all 
y  ^  (  r  ( )  and  y  ^  ( r  .  )  by  this  value,  where  r^  are  the  seven  positions 

at  which  y  ^  and  y  must  be  specified  so  as  to  permit  the  next  step 

of  the  fourth-order  p r ed  ic t o r - c o r r e c t o r  method.  Seven,  rather  than 
four  r  are  required  to  allow  the  automatic  doubling  of  step  size 
when  certain  depths  are  reached. 

For  ease  of  reference,  a  normalization  scheme 
which  is  appropriate  for  the  program  segment  in  Figure  4a,  is 
given  in  Figure  7.  Two  warnings:  (1)  If  only  sparing  use  of 

normalization  is  planned,  or  it  is  to  be  invoked  from  an  IF 
statement,  the  segment  in  Figure  7  will  be  satisfactory;  how- 
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ever,  if  large-scale  use  is  envisioned,  efficiency  requires  the 
inclusion  of  normalization  directly  within  the  coding  in  Figure  4a. 

(2)  Unless  absolutely  necessary,  normalization  should  not  be  in¬ 
cluded  in  these  computations;  it  can  result  in  a  very  significant 
increase  in  computation  time. 

Our  numerical  tests  of  the  overflow  problem  were  performed 
with  the  average  (oceanic)  earth  structure  given  by  Wiggins  (1968). 
The  results  are  given  in  Table  3.  For  IBM  360  equipment,  when  using 

7  0  6  0 

double-precision  computation,  overflow  occurs  when  l_FI  *  10  to  10 

Returning  to  Figure  6,  for  routine  use  of  the  AJP 
formulation,  it  is  important  to  integrate  only  down  to  slightly 
below  the  position  the  lines  in  the  figure,  and  to  thereby 
minimize  computation  time  and  expense.  For  this  purpose  we  have 
devised  empirical  "laws"  for  determining  the  maximum  depth 

H  =  a  -  r0 

to  which  the  integration  must  be  performed  if  the  resulting 
phase  velocity  is  to  be  valid  to  4  significant  figures.  The 
data  for  these  determinations  are  collected  in  Figure  8.  The 
"laws"  specifying  the  number  of  wavelengths  of  structure  to  be 
retained,  if  4  significant  figures  are  desired  in  the  computed 
phase  velocities,  are 

H/A  «  7-c  for  mode  0  (6.01) 

H/A  »9.5-c  for  mode  1  (6.02) 

H/A  -  [11  +  J[_(M-2)]  -  c  for  modes  2-7  (6.03) 

10 


where  A^  is  the  wavelength,  and  M  is  the  mode  number. 
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From  Che  results  shown  in  Figures  6  and  8  we  also 
have  the  maximum  periods  and  values  of  c  which  can  be  determined 
without  taking  the  structure  of  the  core  into  consideration 
(other  than  to  determine  the  values  of  g ( r ) ) .  These  results 
are  shown  in  Figure  9,  along  with  the  corresponding  values  of  the 
minimum  order  number  1 . 

Relative  to  the  computation  of  theoretical  seis¬ 
mograms,  the  combination  of  the  results  in  Figure  9  and  those 
given  by  Schwab  and  Kau-sel  (Section  5,  1976),  indicate  a  potentially 

useful  conclusion:  (1)  When  .  ,  only  the  crust-mantle 

—  min 

system  need  be  used  in  the  computations;  c  can  be  computed  at 
specified,  equally-spaced  frequencies  and  inverse  Fourier  trans¬ 
formation  can  be  used  to  calculate  the  theoretical  seismogram  for 
this  range  of  periods;  and  the  first  term  of  the  asymptotic  ex¬ 
pansions  for  P^m  and  Q^m  can  be  used  (possibly  corrected  by 
automatic  numerical  interpolation  from  the  data  in  Figures  2  and 
3  of  Kausel  and  Schwab  (1976)).  (2)  When  l<l  ,  the  core 

must  be  included  in  the  computations;  u>  should  be  computed  at 
integral  values  of  l_,  and  summation  over  should  replace 
Fourier  synthesis;  and  the  exact,  integral-^  expressions  should 
be  used  to  evaluate  the  associated  Legendre  functions. 
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A  point  of  considerable  interest  to  those  involved 
in  the  actual  computation  of  surface  wave  dispersion,  is  whether 
or  not  a  difficulty  analogous  to  the  Thomson-Haskell  "loss-of- 
precision"  p r ob lem ( S c hwab  and  Knopoff,  1970;  Schwab,  1970)  is 
in^rins. :  to  the  AJP  formulation.  In  computations  based  on 
the  homogeneous - lay e r  approximation,  when  the  original  version 
(Haskell,  1953)  of  the  Thomson-Haskell  formulation  for  Rayleigh 
waves  is  used,  this  problem  can  cause  serious  difficulties  if 
the  computer  is  employed  in  a  low-precision  mode.  To  test  for 
an  analog  to  this  11  loss-of -p r ec is  ion"  problem,  in  our  optimization 
of  the  AJP  formulation,  we  simulated  single-precision  (about  6 
decimal  digits)  computation  by  replacing  DO-loop  160  (Figure  4a) 
in  our  doub le-pr ec is  ion  program,  with  the  program  segment  shown 
in  Figure  10.  The  function  SNGL  accepts  a  double-precision 
argument,  and  returns  the  single-precision  equivalent. 

The  results  of  our  single-precision  tests  are 
similar  to  those  from  the  original  Thomson-Haskell  formulation 
(Figure  2,  Schwab  and  Knopoff,  1970),  and  are  illustrated  in 
Figure  11.  In  the  r-_o  range  shown,  our  results  indicate  that 
there  is  no  problem  with  modes  0-4  in  double-precision  computations 
but  when  computations  are  reduced  to  single  precision,  the  loss-of- 
precision  problem  is  clearly  illustrated.  In  the  latter  case, 
there  is  seen  to  exist  a  minimum  value  of  r_$_,  below  which  we 
cannot  go  and  still  retain  a  given  accuracy  in  the  computed  phase 
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velocities.  Thus  the  AJP  formulation  does  indeed  exhibit  the 
analog  to  the  Thomson-Has ke 11  lo s s -of-pr e cis ion  problem. 

Since  practical  work  with  dispersion  computations 
on  IBM  360  equipment  is  routinely  carried  out  in  double  precision, 
it  is  important  to  estimate  the  numerical  limitations  imposed  by 
the  loss-of -precis  ion  problem  in  this  computational  mode.  The 
results  of  our  tests  at  50  and  25  seconds  are  shown  in  Figures 
12  and  13.  Less  extensive  tests  were  also  carried  out  at  a 
period  of  65  seconds.  At  a  given  period,  the  right-most  point  of 
each  of  the  smoothed  curves  was  used  to  determine  the  maximum 
accuracy  possible  for  each  mode.  This  information  was  then 
collected  in  the  Figures  14  and  15.  Although  the  data  is 
necessarily  sparse,  due  to  the  expense  of  this  type  of  experiment, 
the  results  are  strikingly  clear:  For  a  fixed  period,  as  we  go 
to  higher  and  higher  mode  numbers,  the  attainable  accuracy  in  £ 
becomes  less  and  less;  for  a  fixed  accuracy  in  £,  as  we  go 
to  shorter  and  shorter  periods,  the  maximum  mode  number  that 
can  be  successfully  treated  becomes  smaller  and  smaller. 

In  the  near  future  we  intend  to  present  the  details  of 
our  numerical  analysis  of  the  various  possible  methods  for  dealing 
with  the  loss-of-preciB ion  problem.  See,  for  example,  Wiggins  (1968), 
Nolet  (Appendix  A,  1976),  Nelgauz  and  Skadlnskaya  (1972),  Gilbert 
and  Backus  (1966),  and  Takeuchi  and  Salto  (Section  II. D. 4  ,  1972). 
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TERMINATING  BOUNDARY  CONDITION’S 


35  . 


In  our  present  programming  efforts  we  are  concentrating 
on  the  computation  of  phase  velocities  at  arbitrarily  specified 
(equally-spaced)  frequencies.  To  keep  our  algorithms  as  simple 
as  possible,  and  to  avoid  the  difficulties  involved  in  the 
evaluation  of  spherical  Bessel  functions  of  non-integral  order, 
we  have  chosen  to  terminate  our  integrations  at  depth  with  free, 
or  rigid,  terminating  boundary  conditions.  For  completeness,  we 
also  include  th ■■  technique  for:  (1)  terminating  the  integration 

within  the  mantle  by  applying  terminating  boundary  conditions  for 
a  gravitating,  homogeneous,  solid  sphere  below  rc  and,  (2)  termi_ 
nating  at  the  mantle-core  boundary  by  applying  the  conditions  for  a  homo¬ 
geneous  liquid  sphere  below  r 

0 


In  the  former  case »  just  above  r  we  have 

o 


Y  +  (r(|  )  -  yl(a)X1(r0)  +  y3(b)X2(r0)  +  y5(a)X3(r0)  (7.01) 


where 


a  for  a  continental  structure 


b  - 


(7.02) 


rj  for  an  oceanic  structure 
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For  the  homogeneous,  gravitating,  solid  sphere  below  r^  ,  there 

are  three  classes  of  solutions:  Yi(r),  Y2  (r) ,  Y3  (  r  )  ;  thus  just 

be  low  r,  , 

6 

Y  ( r  q )  =  D  Y  1  ( r  a )  +  £  Y  2  (  r  q  )  +  M  Y  3  (  r  g  )  ,  (7-03) 


where  13,  ji,  and  M  are  undetermined  constants.  Applying  the  boundary 
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or 

NW  -  0  (7.05) 

and  the  dispersion  function  has  the  form 


F  ( co ,  c  )  =  det(N)  . 


(7.06) 
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The  components  of  Y^(r)  are  given  in  convenient  form  by  Takeuchi 

and  Saito  (1972):  Yj(r)  bvtheir  equations  (98)iuith  the  negative 

sign  in  (99);  Y,(r)  by  (98),  with  the  positive  sign  in  (99);  and  Y^r) 

by  their  equations  (100).  Note  that  their  definition  of  y  differs 

6 

slightly  from  that  used  here. 

When  the  structure  used  to  form  the  dispersion  function  is 
terminated  at  the  mantle-core  boundary  oy  the  conditions  for  a 
homogeneous  liquid  sphere  below  r^  ,  above  r^  we  have  (7.01);  below, 

Y_  ( r  <j )  =  P  Y  !  (r  o )  +  Q  V3(r0)  ,  (7.07) 

where  £  and  (J  are  undetermined  constants, 
and  Y  ^  in  (7.07)  have  the  form 

Y.(r0) 

Again,  Takeuchi  and  Saito  (1972)  give  the  form  of  these  vector 
components  for  the  gravitating,  homogeneous,  liquid  sphere.  From 
the  conditions  of  continuity  of  yj,  yz,  V5,  yg  at  r0 , 
vanishing  of  yu(r0),  we  have 


and  the 
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(7.09) 


(7.10) 

(7.11) 

To  obtain  the  group  velocity  we  still  employ  (3.03)  and  (3.04). 

The  forms  of  and  £  which  result  from  (7.06),  or  (7.11),  can  be  obtained  by 
analogy  with  the  way  in  which  (3.05)  is  obtained  from  (2.10).  In  the  present 
case,  however,  the  analog  of  (3.05)  will  comprise  the  sum  of  six,  sixth-order 
determinants  when  (7.06)  is  used  to  form  the  dispersion  function,  and  the  sum 
of  five,  fifth-order  determinants  when  (7.11)  is  used. 

As  would  be  expected,  our  numerical  tests  of  a 
terminating  solid  sphere  show  that,  to  obtain  a  given  accuracy 
in  c  ,  less  structure  must  be  retained  in  this  case  than  when 
using  terminating  rigid  or  free  boundaries  within  the  solid  mantle. 
A  complete  set  of  tests,  comparable  to  those  in  Figure  12,  was 
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and  the  dispersion  function  takes  the  form 

F(u),c)  =  det(R) 
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performed  for  a  period  of  50  seconds;  Che  results  showed  that 
the  increase  in  the  maximum  r_^  values  was  surprisingly  independ¬ 
ent  of  mode  number  and  a  :  (300+25)  km,  or  in  most  cases, 

about  24  fewer  integration  steps  when  terminating  with  a  solid 
sphere.  Thus  the  use  of  the  "correct"  terminating  boundary 
condition  appears  to  be  important  only  for  the  lowest  (radial) 
mode  numbers.  Our  tests  with  mode  7  at  25  seconds,  show  the 
increase  in  r_j_  to  be  about  (110+10)  km.  Thus,  from  our  limited 
number  of  tests,  it  appears  that  g r  p  /  r  p  is  roughtly  constant  with 
about  the  value  0.14+0.02  ;  where  _r_j_  is  maximum  value  required 
by  the  structural  limitation  when  the  condition  of  a  rigid 
terminating  boundary  is  employed,  and  6r  0  is  the  increase  possible 
in  ££  when  one  then  employs  Che  condition  of  a  terminating  solid 
sphere . 

Tests  of  the  los s-o f -p rec i s ion  problem  were  also  per¬ 
formed  with  the  "correct"  boundary  condition  at  depth.  Again,  a  com¬ 
plete  set  of  tests  was  carried  out  at  50  seconds,  and  mode  7  was 
tested  at  a  period  of  25  seconds.  The  right-hand  extremes  of  the 
analogs  of  the  smooth  curves  in  Figures  12  and  13,  occurred  at  the 
same  depths  in  these  new  tests;  thus,  as  a  result  of  the  increased 
r o  values  of  the  upper  portions  of  the  curves  when  solid-sphere 
termination  is  used,  the  maximum  accuracy  for  any  given  mode  is 
significantly  improved  by  using  this  type  of  boundary  at  depth. 
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8 .  CONCLUSIONS 

An  analysis  of  direct,  Rayleigh-wave  dispersion  com¬ 
putations  on  a  spherical,  gravitating  earth  has  been  performed 
using  the  A1 terman-Jarosch-Pekeris  (1959)  formulation. 

1.  No  difficulty  was  encountered  when  we  reversed  the 
usual  procedure,  for  practical  purposes,  and  computed  phase 
velocities  (or  polar  order  numbers)  at  specified  periods. 

2.  Integration  from  the  free  surface  downward,  again 
reversing  the  "standard"  procedure,  resulted  in  no  unexpected 
difficulties.  In  fact,  this  procedure  much  simplified  the 
specification  of  the  algorithm  for  integrating  the  system  of 
differential  equations  to  obtain  phase-velocity  dispersion. 

This  procedure  makes  the  generalization  from  the  algorithm  for 
continental,  to  oceanic  structures  relatively  trivial;  also, 
it  makes  it  possible  to  develop  direct  algorithms  for  obtaining 
group  velocities  for  the  two  types  of  structures.  The  usual 
variational  techniques  are  not  required  for  this  latter  purpose. 

3.  Our  optimization  of  the  AJP  formulation  is  based  on 
removing  all  function  evaluations  from  the  innermost,  integration 
(over  £)  loops  of  the  program.  In  fact,  most  of  the  evaluation 
procedure  for  a ^  ( r ^ )  can  even  be  removed  from  the  phase  velocity 
and  frequency  loops.  This  optimization  of  the  AJP  formulation 
yields  a  characteristic  time  of  336  use conds / in t egr a t ion  step/ 
iteration  for  spherical,  gravitating  structures,  and  a  time  of 
151  useconds / in t eg ra t ion  s t ep / i t era t ion  for  spherical. 


non-gravitating  structures.  These  characteristic  times,  for 
the  AJP  d i r e c t - i n t egr a t ion  procedure,  provide  a  basis  of  com¬ 
parison  with  the  homogeneous- lay e r  approximation  for  a  non¬ 
gravitating  flat  structure,  which  has  a  characteristic  time 
of  110  p se c / lay e r/ i t e ra t ion  (Knopoff's  method,  Schwab  and 
Knopoff,  19  7  2).  The  lowe  r  bounds  of  the  characteristic 
times,  for  our  final  optimizations  of  the  AJP  formulation,  are 
266  psec/step/iteration  for  the  gravitating  case  and 
143  psec/step/iteration  when  gravity  is  not  included.  All  of 
the  above  times  apply  to  the  IBM  360/91  computer  at  UCLA. 

4.  Our  results  h ere , comb i ned  with  those  of  Schwab 

and  Knopoff  (1972),  indicate  that  an  integration  "step"  (in 

the  AJP  procedure)  can  be  considered  nearly  equivalent  to  a  "layer" 

in  computations  based  on  the  homogeneous- 1  ay e r  approximation. 

Also,  to  the  accuracy  possible  in  this  type  of  comparison,  the 
"iterations"  required  in  the  two  techniques  (see  Schwab  and 
Knopoff  (1972)  for  details)  can  be  considered  equivalent.  Thus 
the  relative  efficiencies  of  the  two  types  of  Rayleigh-wave  dis¬ 
persion  computations  can  be  evaluated  by  simple  comparison  of  the 
above  characteristic  times.  The  fact  that  approximately  the  same 
number  and  sizes  of  "steps"  must  be  used  in  the  di re c t- i n t egr a t ion 
procedure,  as  number  and  sizes  of  "layers"  in  the  homogeneous- 
layer  approximation,  means  that  the  usual  assumption  that  the 
former  method  does  a  better  job  of  treating  continuous  parameter- 
depth  distributions,  appears  to  be  invalid. 

5.  The  overflow  problem  in  the  AJP  formulation  can  be 
controlled  by  simple  normalization.  Program  segments  are  given 
which  describe  the  procedure  explicitly. 


6.  A  loss-o f -p r e c is  ion  problem  appears  to  be  intrinsic 
to  the  AJP  formulation.  Results  of  this  problem:  For  a  fixed 
period,  as  mode  number  increases,  the  attainable  accuracy  in 
the  phase  velocity  decreases;  for  a  fixed  accuracy  in  the  phase 
velocity,  as  period  decreases  the  maximum  mode  number  that  can  be 
treated  successfully  decreases. 
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Table  1.  Constants  for  integration  through  successive  step-size  regions 
illustrated  in  program  segment  given  in  Figure  4.  Constants 

correspond  to  4  significant  figures  in  computed  phase  velocity. 

I  N1(I)  N2(I)  Integration 

step  size 
(km) 

1  5  11  -1.5625 

2  12  16  -3.1250 

3  17  21  -6.2500 

4  22  *  -12.5000 

*  N2 ( 4 )  is  specified  so  as  to  allow  integration  to  proceed  to  the 
deepest  point  within  the  solid  mantle,  while  maintaining  a  step 
size  of  -12.5  km.  NEND  is  determined,  at  each  period,  by  the  input 
value  of  r„;  it  must  satisfy  NEND  ^  N2(4). 


able  2.  Relative  speeds,  for  surface-wave  dispersion  computations,  for  several  computers. 
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Table 


3.  Results  of  numerical  tests  of  the  overflow 

problem  when  normalization  is  not  included  in  our 
optimization  of  the  basic  AJP  formulation.  An 
average  (oceanic)  earth  structure  (Wiggins,  1968), 
and  a  period  of  50  seconds,  were  used  in  the  tests. 
The  value  of  r0  is  the  maximum  at  which  overflow 
occurs;  H/X  is  the  corresponding  number  of  wave¬ 
lengths  of  structure,  from  the  surface  of  the  earth 
down  to  r=*  r  o  .  These,  and  larger  values  of  H/X  ,  yield 
overflow. 
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Figure 


Schematic  representation  of  optimized  scheme  for 
evaluating  the  matrix  elements  a^  ^ ( r in  the 
treatment  of  the  solid  sedimentary  layers,  the 
subsedimentary  crustal  layers,  and  the  mantle. 

The  quantities  and  p  are  Lame's  constants,  G 
is  the  gravitational  constant,  and  is  the  number 
of  depths  at  which  a  (r^)  oust  be  evaluated. 


Schematic  representation  of  optimized  scheme  for 

evaluating  the  matrix  elements  b,  (t  )  in  the 

_ _ _ 

t rea tmen t  of  the  homogeneous  oceanic  (liquid)  layer. 
The  quantity  a  Is  the  compressional-wave  velocity, 

£  is  the  density,  g (r )  is  acceleration  due  to  gravity, 
and  M  is  the  number  of  depths  at  which  b^  (r^)  must 
be  evaluated. 


3.  (a)  FORTRAN  IV  program  segment  for  the  basic  matrix 

multiplication  for  our  optimization  of  the  AJP 
formulation  for  solid  layers;  (b)  symbolic  repre¬ 
sentation  of  (a),  which  is  used  in  Figure  4;  and 
(c)  definition  of  one-dimensional  array  used  in  (a) 
to  represent  6x6  matrix  in  (2.01).  The  integer  IPT 
i s  the  index  specifying  the  value  of  r_,  and  A1 
through  A36  are  dimensioned  to  300. 

4a.  FORTRAN  IV  program  segment  in  which  the  predictor- 

corrector  portion  of  the  integration  from  belov  the  Moho 
to  ro  is  handled;  most  of  the  computation  time  is  spent 
in  this  segment.  The  boxed  segments  refer  to  the  basic 

S  TT*  -...Is-*..  IT--.!  Ml,.  ~  t  -» 


Figure  4b.  FORTRAN  IV  program  segment  demonstrating  sub¬ 
scripting  and  storage  improvements,  relative 
to  the  segment  in  Figure  4a,  that  are  required 
to  optimize  computation  time  on  an  IBM  360 
compu  ter . 

Figure  4c.  FORTRAN  IV  program  segment  illustrating  our 

final  optimization  for  the  non-gravitating  case. 


Figure  5  . 


Oceanic  and  continental  (shield)  models  used  in 
program  testing. 


Figure  6. 


Figure  7  . 


Figure  8. 


Figure  9  . 


Values  of  j^_,  the  depth  at  which  integration  is 
terminated,  which  yield  4-s ignif ican t- f  igur e 
accuracy  in  the  computed  values  of  c_  with  the 
optimized  version  of  the  basic  AJP  formulation. 

At  each  period,  4-figure  accuracy  is  attained 
only  if  r „  is  specified  to  be  smaller  than  the 
indicated  curve  . 

Normalization  scheme  appropriate  to  program  segment 
in  Figure  4a.  The  procedure  should  be  included 
between  statement  numbers  160  and  170  in  Figure  4a. 
See  text  for  warnings  concerning  loss  of  efficiency 
when  normalization  is  employed. 


Results  of  tests  for  determining  general,  multimode 
"laws"  for  specifying  the  required  values  of  r0  to 
use  when  computing  phase-velocity  dispersion  for 
mantle  Rayleigh  waves,  if  c  is  desired  to  4-figure 
accuracy. 


Maximum  periods,  maximum  phase  velocities,  and 
minimum  order  numbers,  f  ,  that  can  be  used  if 
c  is  desired  to  4-figure  accuracy,  when  the  inte¬ 
gration  is  limited  to  the  mantle  and  the  core  is 
excluded  fr  a  the  computations  (other  than  for  use 
in  the  determination  of  g ( r )  )  . 


Figure  iO. 


FORTRAN  IV  Program  segment  used  to  simulate 


s ingle-p rec is  ion  computations  when  using  our 
double-precision  optimization  on  the  basic  AJP 
formulation.  This  program  segment  is  used  to 
replace  DO-loop  160  in  Figure  4a. 

Figure  11.  Test  of  the  effect  of  reducing  IBM  360  com¬ 
putations  from  double  (about  16  decimal  digits) 
to  single  precision  (about  6  decimal  digits), 
while  keeping  period  fixed  at  50  seconds.  For 
a  given  mode,  in  order  to  obtain  a_  significant 
figures  in  the  computed  phase  velocity,  r_o_  must 
not  exceed  the  value  given  by  the  upper  portion 
of  the  dashed  line  (structural  limitation)  nor 
fall  below  the  lower  portion  ( loss-of -pr e c is  ion 
limitation)  . 

Figure  12.  IBM  360  double-precision  tests  of  loss-of -p rec is  ion 
problem  at  a  period  of  50  seconds.  At  left,  raw 
results  of  structure  reduction  experiments;  at  right, 
smoothed  curves  for  each  mode.  Latter  curves  are 
drawn  such  that  all  data  points,  for  a  given  mode, 
fall  to  the  right  of  corresponding  curve. 

Figure  13.  IBM  360  d oub 1 e- p r e c i s ion  tests  of  los s-of -p r ec is  ion 


problem  at  a  period  of  25  seconds. 


Figure  14. 


for  which 


Maximum  possible  mode  number,  n  , 

max 

o  significant  figures  can  be  obtained  vith  our 
optimization  of  the  AJP  formulation.  This 
limitation  is  due  to  the  loss-of -p r ec i s ion 
problem. 

Figure  15.  Relationship  between  maximum  possible  mode  number 
n  ,  and  minimum  period,  for  several  values  of  a 
This  limitation  of  our  optimization  of  the  AJP 
formulation  is  due  to  the  loss-of-precision  probl 


Elements 

to  be  evaluated 

external 

to  both  to  and  c  loops 

a  1  1 

a  1  2 

a  3  3 

a  4  5 

a*  1 

a2  2 

a  3  4 

a  2  6 

as  1 

a  m  2 

a  4  4 

a  6  6 

Auxiliary  quantities  to  be  evaluated  external  to  both  w  and 
c  loops : 

2 

d  2 1  =  4y ( 3X+2y ) aj 2/r 
2 

d„  3=-2y/r 
dj  3  =  Xai 2/r 
d6  3  =  -4ttGp 
d65=l/r* 

2 

e<.  3=4y  (X+y)  ai  2/r 
(a  loop 

2 

OMEGSQ=w 
DO  10  1=1, N 
TEMP=-a2  6 (I) *0MEGSQ 
a2 i  (I) =d2  j  ( I ) -TEMP 
10  f 4  3  (I) =d„ 3 (I) -TEMP 
c  loop 

ORDER=MH+l) 

DO  20  1=1, N 
a i 3 ( I ) =d i 3 (I) xORDER 
a2 3 (I)=aH i (I) xORDER 
a6  3 (I)=d6  j (I) xORDER 
a2  * ( I ) =a 3  j (I) xORDER 
a$  5 ( I ) =ds  s(I) xORDER 


20  a«  3  (I)  =  ft,  3  (I)  +e*  3  (I)  xORDER 


Evaluate  b66  external  to  both  id  and  c  loops 
quantities  to  be  evaluated  external  to  both 
2 

XLAINV=l/a  p  p6s=P&2P  P21- 


2 

p, 2=-l/pr 


2 

Pt  2=4uG/r 
2 


Pi s=-l/r 


P6i="P65g(r)  s21- 

2 

p22=-g(r)/r 
p  2  s  =P  2  2  P 


-  o  loop 

2 

RHMOSQ=-p(i) 

2 

OMSQIN=l/a) 

DO  10  1=1, M 
h2 1  ( I ) =P  2 1 (I) xOMSQIN 
qj  1  (I)=S2 1  (I) +RHMOSQ 
h6 1  ( X ) =P  6 1  ( I ) xOMSQIN 
h:  2  (I)=Pi 2  (I) xOMSQIN 
hJ2 ( I )  =p  2  2 (I) xOMSQIN 
h 6  2  (I)=P$  2 (X) xOMSQIN 
hj  5 ( I ) =pi 5 (I) xOMSQIN 
h25(I)=p25(I)xOMSQIN 
10  h6  5 (I) =P6  s(D  xOMSQIN 
p-  c  loop 

ORDER=i.  U+l) 

DO  20  1=1, M 

b2  2 (I) =h2 2 (I) xORDER 

bi  ,  ( I ) =b$  6 (I) -b2  2 (I) 

b  2 1  ( I ) =h  21(1) xORDER+q  2 J (D 

b  6 1  ( I ) =h  6 1 ( I ) x ORDER 

bi 2 (I) =hi 2 (I) xORDER+XLAINV 

bg2 (I)=b&2(iJ xORDER 

b, s (I)=h, 5 (I) xORDER 

b?  s ( I ) =h2  s ( I ) xORDER 


Auxiliary 
ui  and  c  loops 

P2  sg (r) 

4 ti g  (r)  /r 


U/ 


(b) 


(c) 


M- 


YBAR(l)  : 
YBAR(2) = 

1 

YBAR ( 3) ; 
YBAR ( 4 )  = 

1 

YBAR ( 5) = 
YBAR (6) = 


=A1  ( IPT) *Y(1) +A7 (IPT) *Y(2) +A13 (IPT) *Y (3) 

;A2 (IPT) *Y (1 ) +A8 (IPT) *Y (2) +A14 (IPT) *Y(3)  + 
A20 (IPT) *Y (4) +A32 (IPT) *Y(6) 

:A15 (IPT) *(  Y ( 3 ) -Y ( 1) ) +A21 (IPT) *Y(4) 

:A4 (IPT) *Y ( 1 ) +A10 ( IPT) *  Y ( 2 ) +A16 (IPT) *Y(3)+ 
A22 (IPT) *Y (4) +A2  8 (IPT) *Y  (5) 

;A5 ( IPT) *Y ( 1) +Y  ( 6 ) 

:A18 ( IPT) *Y { 3) +A30 (IPT) *Y ( 5) +A36 ( IPT) *Y(6) 


[ 

YBAR, 

IPT 

] 

A1 

A7 

A13 

0 

0 

0 

A2 

A8 

A14 

A20 

0 

A32 

-A15 

0 

A15 

A21 

0 

0 

A4 

A10 

A16 

A22 

A28 

0 

A5 

0 

0 

0 

0 

1 

0 

0 

A18 

0 

A30 

A36 

C  BEGIN  APPLICATION  OF  PREDICTOR-CORRECTOR  METHOD. 

DO  110  1=1,6 
110  PMNUSC { I ) =0 . OD+OO 

HH=0 . 5D+00  *  1 . 5625D+00 

C  LOOP  OVER  REGIONS  WITH  DIFFERENT  STEP  SIZES. 

DO  180  I REG= 1 , NUMREG 

HH=HH+HH 

NSTART=N1 (IREG) 

NSTOP=N2 (IREG) 

NTEMP=5 

IT=0 

C  LOOP  OVER  DEPTH  IN  CURRENT  STEP-SIZE  REGION. 

DO  170  N=NSTART , NSTOP 
IF ( IT. EQ. 4)  GO  TO  115 
IT=NTEMP-4 
ITP1=IT+1 
ITP3=IT+3 
ITP8=IT+8 
ITP9=IT+9 
ITP10=IT+10 
115  DO  120  1=1,6 
C  SET  PREDICTOR  P(I). 

P (I) =B ( IT , I) +COEFF1* ( 2 . OD+00 * (B (ITP10 , 1 ) +B(ITP8,I)  ) 
-B ( ITP9 , I ) ) 

C  SET  MODIFIED  PREDICTOR  XM(I). 

120  XM ( I) =P ( I ) - . 92561983471074  38D+00*PMNUSC ( I ) 

IPT=IPT+1 


||  XMBAR,  IPT  || 

DO  130  1=1,6 
C  SET  CORRECTOR  C(I) 

C (I) =. 125D+00* (9 . 0D+00*B ( ITP3 , I ) -B ( ITP1 , I) 
+COEFF2  * ( XMBAR ( I ) +2 . OD+00  *B ( ITP10 , 1) 
-B ( ITP9 , I) ) ) 

PMNUSC ( I ) =P ( I ) -C (I) 

C  SET  SOLUTION  VECTOR  AT  NTH  DEPTH. 

130  Y ( I ) =C ( I ) + . 0  7  4  38016  52  89  2  5620D+00 
*PMNUSC ( I ) 

IF(N.EQ.NEND)  RETURN 

C  SET  DERIVATIVE  OF  SOLUTION  VECTOR  AT  NTH  DEPTH. 


||  YBAR,  JPT  || 

IF (NTEMP.GT. 7)  GO  TO  150 
N'j  MPP7=NTEMP+7 
DO  140  1=1,6 
B ( NTEMP , I ) =Y ( I ) 

140  B (NTMPP7 , I)=YB^R(I) 
NTEMP=NTEMP+1 
GO  TO  170 
150  DO  160  1=1,6 
B ( 1 , 1) =B (2,1) 

B ( 2 , I ) =B (3,1) 

B ( 3 , I ) =B (4,1) 

B(4 , I) =B (5, I) 


160 

170 


B (5 , 1) =B (6,1) 

B (6 , 1) =B ( 7 , 1 
B  ( 7 , 1 ) =Y ( I ) 

B ( 8 , 7) =B (9 , X) 

B (9 , I ) =B (10,1) 

B ( 10 , I ) =B (11,1) 

B ( 11 , I) =B (12,1) 

B (12 , I) =B (13,1) 

B ( 13 , I) =B (14,1) 

B ( 14 , I ) =YBAR ( I ) 

CONTINUE 

C  RESET  STORED  VALUES  OF  COEFFICIENTS  IN  PREPARATION  FOR 
C  DOUBLED  STEP  SIZE. 

COEFFl=COEFFl+COEFFl 

COEFF2=COEFF2+COEFF2 

COEFF6=CO£FF6+COEFF6 

C  RESET  STORED  VALUES  OF  Y ( I ) , YBAR ( I ) , and  PMNUSC(I)  IN 
C  PREPARATION  FOR  DOUBLED  STEP  SIZE. 

DO  180  1=1,6 

PMNUSC (I) =8. 962962962962963 D+ 00* (Y(I) 

-B (1,1)  ) -COEFF6  * (YBAR (I) +B(8,I) 

+  3 . 0D+00  * (B ( 12 , 1) +B ( 10 , 1 ) ) ) 

B (2 , I) =B (3,1) 

B  (  3 , 1 ) =B (5,1) 

B ( 4 , 1) =B (7,1) 

B ( 9 , 1) =B (10,1) 

B  ( 10 , 1 ) =B (12,1) 

180  B  ( 11 , 1 ) =B (14,1) 


DIMENSION  B(6,14)  ,Y(6)  ,YBAR(6)  ,XM(6)  ,XMBAR(6)  ,P(6) ,C(6) ,PMNUSC(6)  , 
1  A(2090) 

EQUIVALENCE  (Y(1),Y1) 

EQUIVALENCE  (XM(1) ,XM1) 

EQUIVALENCE  (XMBAR(l) .XMBARl) 

IPT=-18 


C  LOOP  OVER  DEPTH  IN  CURRENT  STEP-SIZE  REGION. 

DO  170  N=NSTART ,NSTOP 
IF ( IT . EQ . 4)  GO  TO  115 
IT=NTEMP-4 
115  DO  120  1=1,6 
C  SET  PREDICTOR  P (I)  . 

P(I)=B(I, IT)+C0EFF1* ( 2 . OD+OO* (B(I , IT+10)+B (I , 1T+8) ) -B(I , IT+9) ) 

C  SET  MODIFIED  PREDICTOR  XM(I). 

120  XM(I)=P(I)-.9256198347107438D+00*PMNUSC(I) 

IPT=IPT+19 

XMBAR(1)=A(IPT)*XM(1)+A(IPT+1)*XM(2)+A(IPT+2)*XM(3) 

XMBAR  ( 2) =A (IPT+3) *XM(l)+A(IPT+4) *XM(2)+A(IPT+5) *XM(3) 

1  +A( IPT+6) *XM(4)+A(IPT+7) *XM(6) 

XMBAR(3) =A(IPT+8) * (XM (3) -XM(1) )+A(IPT+9) *XM(4) 
XMBAR(4)=A(IPT+10)*XM(1)+A(IPT+11)*XM(2)+A(IPT+12)*XM(3) 

1  +A(IPT+13)*XM(4)+A(IPT+14)*XM(5) 

XMBAR(5)=A(IPT+15)*XM(1)+XM(6) 

XMBAR (6) =A (IPT+16) *XM(3)+A(IPT+17) *XM(5)+A(IPT+18) *XM(6) 

DO  130  1=1,6 
C  SET  CORRECTOR  C(I)  . 

C( I) = • 125D+00* (9 . 0D+00*B (I , IT+3) -B (I , IT+1)+C0EFF2* (XMBAR (I) +2 . 0D+0 
1  *B(I , IT+10) -B (I , IT+9) ) ) 

PMNUSC (I) =P (I)-C (I) 

C  SET  SOLUTION  VECTOR  AT  NTH  DEPTH. 

130  Y ( I) =C (I )+. 07438016528925620D+00*PMNUSC (I) 

IF(N. EQ.NEND)  RETURN 

C  SET  DERIVATIVE  OF  SOLUTION  VECTOR  AT  NTH  DEPTH. 

YBAR(l) =A ( IPT) *Y ( 1)+A( IPT+1) *Y (2)+A (IPT+2) *Y (3) 

YBAR(2) =A( IPT+3) *Y (1) +A( IPT+4) *Y (2)+A(IPT+5) *Y ( 3) 

1  +A( IPT+6) *Y (4)+A(IPT+7) *Y (6) 

YBAR ( 3) =A( IPT+8) * (Y (3) -Y ( 1) )+A(IPT+9) *Y (4) 
YBAR(4)=A(1PT+10)*Y(1)+A(IPT+11)*Y(2)+A(IPT+12)*Y(3) 

1  +A( IPT+13) *Y(4)+A( IPT+14) *Y (5) 

YBAR (5) =A(IPT+15) *Y (1)+Y (6) 

YBAR ( 6) =A( IPT+16) *Y ( 3)+A(IPT+17) *Y (5)+A(IPT+18) *Y (6) 

IF(NTEMP.GT. 7)  GO  TO  150 
DO  140  1=1,6 
B ( I , NTEMP ) =Y ( I ) 

140  B(I .NTEMP+7) =YBAR(I) 

NTEMP=NTEMP +1 
CO  TO  170 


C  RESET  STORED  VALUES  OF  Y  AND  YBAR  IN  PREPARATION  FOR  NEXT  INTEGRATION 
C  STEP. 

150  DO  16C  1=1,6 
B(I,1)»B(I,2) 

B(I , 2)=B(I , 3) 

B(I,3)**B(I,4) 

B(I,4)=B(I,5) 

B(I ,5)=B(I ,6) 

B(I , 6)=B(I , 7) 

B(I,7)=Y(I) 

B(I,8)=B(I,9) 

B(I,9)=B(I, 10) 

B (I , 10)=B(I ,11) 

B(I , 11)=B(I , 12) 

B(I , 12) =B (I , 13) 

B(I , 13)=B(I , 14) 

160  B (I , 14) =YBAR(I) 

170  CONTINUE 


DIMENSION  B(4,]4) ,D(4,14) ,X{4) ,Y(4) ,XBAR(4) ,YBAR(4) ,XM(4)  ,YM(4)  , 

1  XMBAR(4) ,YMBAR(4) ,P(4) ,Q(4) ,C(4) ,F(4) ,PMNUSC(4) , 

2  QMNUSF(4) ,A(1430) 

EQUIVALENCE  (X(1),X1) 

EQUIVALENCE  (Y(l) ,Y1) 

EQUIVALENCE  (XM(1) ,XM1) 

EQUIVALENCE  (YM(1) ,YM1) 

EQUIVALENCE  ( XMBAR ( 1 ) , XMBAR 1 ) 

EQUIVALENCE  (YMBAR(l) YMBAR1) 

IPT=-12 


C  LOOP  OVER  DEPTH  IN  CURRENT-STEP-SIZE  REGION. 

DO  170  N=NSTART ,NSTOP 
I F ( IT. EQ . 4)  GO  TO  115 
IT=NTEMP-4 
115  DO  120  1-1,4 

C  SET  PREDICTORS  P(I)  ,  Q (I)  . 

P(I)=B(1, IT)+C0EFF1*(2 . OD+OO* (B(I , IT+10)+B(I , IT+8) )-B(I , IT+9) ) 
Q(I)=D(I,IT)+COEFFl*(2.0D+00*(D(I,IT+10)+D(I,IT+8))-D(I,IT+9)) 

C  SET  MODIFIED  PREDICTORS  XM(I)  ,  YM(I)  . 

XM(I)=P(I)-.9256198347107438D+00*PMNUSC(I) 

120  YM(I)=Q(I)-. 9256198347107438D+00*QMNUSF(I) 

IPT=IPT+13 

XMBAR (1)  =A(1PT)  *XM(1)+A(IPT  t-1)  *XM(2)+A(IPT+2)  *XM  (3) 

YMBAR (1)=A( IPT) *YM (l)+A(IPT+l)*YM(2)+A(IPT+2) *YM(3) 

XMBAR ( 2) =A ( I PT+3) *XM ( 1 ) +A ( IPT+4) *XM  ( 2) +A ( IPT+5) *XM ( 3) 

1  +A( IPT+6) *XM(4) 

YMBAR(2)=A(IPT+3)*YM(1)+A(IPT+4)*YM(2)+A(IPT+5)*YM(3) 

1  +A(IPT+6)*YM(4) 

XMBAR (3) =A( IPT+7) * (XM(3) -XM(1) )+A(IFT+8) *XM(4) 

YMBAR(3)=A(IPT+7)  *  (Y'M{3)  -YM(1)  )+A  (IPT+8)  *YM(4) 

XMBAR(4) =A (IPT+9) *XM( 1)+A(IPT+10) *XM (2)+A(lPT+ll) *XM(3) 

1  +A(IPT+12)*XM(4) 

YMBAR (4) =A( IPT+9) *YM(1)+A(IPT+10) *YM(2)+A(IPT+11) *YM(3) 

1  +A(IPT+12)*YM(4) 

DO  130  1=1,4 

C  SET  CORRECTORS  C(I)  ,  F(I)  . 

C(I)  = .  125D+00* (9 . 0D+00*B( I , IT+3)-B (I , IT+l)+COEFF2* (XMBAR (I)+2 .OD+O 
1  *B(I , IT+10) -B(I , IT+9) ) ) 

F(I)  = . 125D+00* (9 . 0D+00*D(I , IT+3)-D(I , IT+l)+COEFF2*(YMBAR(I)+2  .OD+O 
1  *D(I , IT+10)-D(I, IT+9)) ) 

PMNUSC(I)=P(I)-C(I) 

QMN'USF(I)  =Q(I)  -F  (I) 

C  SET  SOLUTION  VECTORS  AT  NTH  DEPTH. 

X( I) =C(I)+. 074 380165289 25620D+00*PMNUSC(I) 

130  Y(I) =F(I)+.07438016528925620D+00*QMNUSF(I) 

IF(N.EQ.NEND)  RETURN 


I 


C  SET  DERIVATIVES  OF  SOLUTION  VECTORS  AT  NTH  DEPTH. 

XBAR(l) =A( IPT) *X( 1)+A(IPT+1) *X(2)+A(IPT+2) *X (3) 

YBAR(l) =A(IPT) *Y (1)+A (IPT+1) *Y(2)+A(IPT+2) *Y (3) 

XBAR(2) =A( IPT+3) *X ( l)+A(IPT+4) *X(2)+A(IPT+5)*X(3) 

1  +A(IFTP6)*X(4) 

YEAR (2) -A (IPT+3) *Y (1)+A(IPT+4)*Y (2)+A(IPT+5)*Y (3) 

1  +A(IPT+6) *Y (4) 

XBAR (3) =A(IPT+7) * (X(3) -X(l) )+A(IPT+8) *X(4) 

YBAR (3) =A( IPT+7) * (Y (3) -Y (1) )+A{IPT+8) *Y (4) 
XBAR(4)=A(IPT+9)*X(1)+A(IPT+10)*X(2)+A(IPT+11)*X(3) 

1  +A( IPT+12) *X (4) 

YBAR(4) =A( IPT+9 ) *Y (1)+A(IPT+10) *Y (2)+A(IPT+ll) *Y(3) 

1  +A(IPT+12)*Y (4) 

IF(NTEMP.GT.7)  GO  TO  150 
DO  140  1=1,4 
B(I,NTEMP)=X(I) 

D(I,NTEMP)=Y(I) 

B(I,NTEMP+7)=XBAR(I) 

140  D(I,NTEMP+7)=YBAR(I) 

NTEKP=NTEMP+1 
GO  TO  170 

C  RESET  STORED  VALUES  OF  X  ,  Y  AND  XBAR  ,  YBAR  IN  PREPARATION  FOR  NEXT 

C  INTEGRATION  STEP. 

150  DO  160  1=1,4 
B(1,1)=B(1,2) 

B(I,2)=B(I,3) 

B(I,3)=B(I,4) 

B(I,4)=B(1,5) 

B(I ,5)=B(1 ,6) 

B ( I ,6)=B(I ,7) 

B(1 , 7)=X(I) 

B(I,8)=B(I,9) 

B(I,9)=B(I, 10) 

B(I , 10)=B(I ,11) 

B(I , 11)=B(I , 12) 

B(1 ,12)=B(1 , 13) 

B(I , 13) =B (I , 14) 

B (1 , 14)  =XBAR(I) 

D(I,1)=D(I,2) 

D(I,2)=D(I,3) 

D(I,3)=D(I,4) 

D(I,4)=D(X,5) 

D(I,5)=D(I,6) 

D(I,6)=D(I,7) 

D(I,7)=Y(I) 

D( 1 , 8) =D(1 ,9) 

D(I , 9) =D(I , 10) 

D(1 , 10)  =D(1 , 11) 

D (I ,11) =3(1,12) 

D ( I , 12) =D(I , 1 3) 

D(I , 13)=D(1 , 14) 

160  D(I,14)=YBAR(I) 

170  CONTINUE 


C  NORMALIZATION. 

AMXINV=1 . 0D+00/DMAX1 ( DABS (B ( 7 , 1 ) ) , DABS ( B ( 7 , 2 ) ) , 

1  DABS (3(7,3)) , DABS ( B ( 7 , 4 ) ) , 

2  DABS (B (7 , 5) )  ,  DABS (B  (7 , 6)  )  ) 
DO  165  1=1,6 

PMNUSC(I) =PMNUSC (I) *AMX1NV 
Y(I)=Y (I) *AMXINV 
YBAR(I) =YBAR(I) *AMXINV 
DO  165  J=1 , 14 
165  B(J,I)=B(J,I) *AMXINV 


B(1,I)«SNGL(B(2,I)) 
B(2,I)=SNCL(B(3,I)) 
B(3,I)=SNGL(B(4,1)  ) 
B(4,I)=SNGL(B(5,I)) 

B ( 5 , I ) =SNGL ( B  (  6  ,  I)  ) 
B(6,I)=SNGL(B(7,I)) 

B( 7, I) =SNGL(Y (I) ) 
B(8,I)=SNGL(B(9,I)) 
B(9 , I) =SNGL(B(10,  I)  ) 
B(10,I)=SNGL(B(11,I)) 
B(11,I)=SNGL(B(12,I)) 
B(12,I)=SNGL(B(13,1)) 
B(13,I)-SNGL(B(14,I)) 
160  B(14,I)=SNCL(YBAR(I)) 
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RESEARCH  NOTE 

GENERATION  OF  COMPLETE  THEORETICAL  SEISMOGRAMS 
FOR  SH  III 
Enzo  Mantovani 


Summary 

Earlier  efforts  to  generate  the  entire  theoretical 
seismograms,  including  both  body  and  surface  waves  for 
realistic  sources  buried  in  a  radially  heterogeneous 
anelastic,  spherical  Earth,  are  extended  to  include  the 
summation  of  sixteen  modes.  The  comparison  between  a 
real  seismogram  and  theoretical  time  series,  relative  to 
different  attenuation  models  in  the  upper  mantle,  yields 
information  concerning  the  anelasticity  under  the  Pacific 


Ocean . 


Introduction 

Complete  seismograms  which  include  both  body  and  surfg 
ce  waves  for  a  spherical,  anelastic,  radially  heterogeneous 
Earth  have  been  generated  by  simple  inverse  Fourier  tran¬ 
sformation  of  the  propagating  fundamental  and  higher-mode 
surface  waves  (Nakanishi  et  al.,1977,  Mantovani  et  al. ,1977). 
Efficient  computational  algorithms  permit  us  to  deal  with 
highly  realistic  models  of  the  Earth,  whose  radial  hetero¬ 
geneity  is  approximated  by  200  layers,  and  to  employ  2000 
frequency  points  per  radial  mode. 


Theoretical  seismograms 

The  tests  we  report  here  were  carried  out  on  a  CIT-1,1. 
oceanic  structure  (Anderson  and  Toksoz,  1963).  ^he  first 
sixteen  Love  wave,  or  torsional  modes  were  used,  and  for  e- 
ach,  the  dispersion,  attenuation  and  excitation  are  computed 
down  to  a  minimum  period  of  1  second.  The  source  model  we 
have  adopted  is  a  dip-slip  displacement  dislocation  on  a  ve£ 
tlcal  fault  at  a  depth  of  180  km. 

In  Fig.  1  we  illustrate  the  improvement  in  the  details  of 
the  theoretical  seismograms  as  the  number  of  radial  modes 
used  in  the  synthesis  increases.  From  top  to  bottom,  as  the 
number  of  modes  increases, the  amplitude ,  clarity  and  short 
period  content  are  seen  to  increase  for  each  arrival  in  a 
manner  corresponding  to  the  depth  of  penetration  of  the  mo¬ 
des  and  the  depth  at  which  the  arrival  is  produced:  the  first 
six  modes  appear  sufficient  to  reproduce,  almost  completely, 
the  Sa  surface  wavea  train,  that  is,  the  late-arriving,  domJL 
nant  energy  whose  penetration  is  limited  to  approximately 
400km  of  depth  (Kausel  et  al.,1977);  in  the  second  trace, 
synthesized  from  the  first  eleven  modes,  also  the  shallow 
body  wave  phase  SS,  having'  a  group  velocity  of  about  5.7 
km/sec*  appears,  to  be  almost  completely  formed;  in  the  last 
trace  ,  synthesized  from  16  modes,  the  deeper  body  wave  phg 
sea  S  and  SS,  l.e.,  the  first  bursts  on  the  time  series  be¬ 


gin  to  appear  as  two  resolvable  pulses,  but  are  still  devel£ 
ping  injmplitude  and  short-period  content. 


The  core  reflections  ScS  and  sScS,  whose  penetration  is  much 
greater  than  that  of  the  first  16  modes,  at  the  predominant 
periods  comprising  these  pulses,  are  still  unrecognizable 
in  the  seismograms. 

In  Fig.  2  we  show  a  range  of  results  out  to  an  epic.n_ 
tral  distance  of  8000  km;  the  direct  waves  and  those  reflec¬ 
ted  by  the  surface  of  the  Earth  are  recognizable  in  almost 
every  case;  the  waves  reflected  by  the  core  however  are  re¬ 
presented  only  by  their  long  period  components  (  see  traces 
for  distances  of  3000  and  4000  km.).  Relative  to  the  corre¬ 
sponding  figure  in  Part  II  (Mantovani  et  al.,1977),  it  is 
possible  to  note  a  general  improvement  in  the  Impulsive 
shape  and  resolution  of  individual  body  waves  and  is  parti¬ 
cularly  clear  at  short  distance  in  those  traces  for  the 
15-100  WWSSN  instrument. 

The  comparison  between  experimental  and  theoretical  sei 
smograms  for  30-100  WWSSN  instruments  is  given  in  Fig. 3. 

Details  concerning  the  event  and  a  preliminary  compari¬ 
son  between  theoretical  and  experimental  time  series  using 
surface  waves  is  reported  by  Kausel  et  al.  (1977).  Our  pur¬ 
pose  here  is  to  continue  the  discussion,  begun  in  Part  II, on 
the  qualitative  aspects  of  the  theoretical  -  experimental 
comparison  as  more  and  more  higher  modes  are  included  in  the 
synthesis  of  the  theoretical  traces.  We  also  report  on  some 
tests  of  the  sensitivity  of  the  amplitudes  and  arrival  times 
of  the  pulse  bursts  to  the  model  of  the  Earth's  intrin 

sic  attenuation  and  to  the  source  parameters. 

In  Part  II  the  variation  of  relative  amplitudes  of  body 
and  surface  waves  caused  by  changes  in  crustal  and  upper-man 
tie  attenuation  was  noted,  and  it  was  sugg^-ted  that  this 
effect  might  be  a  useful  tool  in  the  investigation  of  the 
Earth's  intrinsic  anelasticity.  We  have  used  two  models  of 
attenuation.  Model  I  is  based  on  model  K1C8  of  Anderson  et 
al.(l965).  This  model  was  obtained  from  measurements  of  the 
decay  of  Rayleigh  and  Love  waves  traveling  globe  circling 
paths;  thus  it  cannot  be  expected  to  be  an  accurate  represen¬ 
tation  of  the  anelasticity  in  a  specific  oceanic  region. 


Model  II  Is  based  on  that  given  by  Mitchell  (1976)  derived 
from  inversion  of  Rayleigh-wave  attenuation  measurements  in 
the  "older  portions  of  the  Pacific  crust  and  upper  mantle". 
Models  I  and  II  are  shown  in  Fig. 4. 

Since  the  average  anelasticity  in  the  first  400  km  is  very 
similar  for  model  I  (  Q  =  0.00825  )  and 

model  II  (  Q  =  0.00774  )  the  Sa  waves 

have  about  the  same  attenuation.  The  ratio  of  amplitudes  of 
Sa/S  is  about  the  same  for  both  models  since  amplitudes  of 
S  are  only  slightly  affected  by  variations  of  anejasticity 
in  a  small  portion  of  their  path.  Since  Sa  waves  sample  the 
uppermost  400  km  of  the  structure  almost  uniformly,  they  are 
more  sensitive  to  the  average  attenuation  in  this  region  than 
to  its  detailed  depth  distribution. 

A  comparison  of  traces  A  and  B  (Fig. 3),  with  the  help  of  TA¬ 
BLE  1,  indicates  that,  for  a  source  at  180  km,  the.  relative 

amplitudes  of  S  and  Sa  are  not  so  different  for  model  I  and 

~  “■  ■■■  • 

II.  The  striking  feature  of  trace  B,  relative  to  A,  is  the 
remarkable  increase  in  the  amplitudes  of  SS  arrivals.  It  is 
reasonable  to  suppose  that  this  effect  is  connected  with  the 
lower  attenuation  of  model  II  between  about  250  and  700  km, 
that  is  in  that  part  of  the  structure  in  which  the  SS  rays 
have  a  elgnif leant  portion  of  their  path. 

A  comparison  of  the  well  pronounced  SS  phase  in  trace  B,  with 
even  the  largest  amplitudes  around  the  same  arrival  in  the 
experimental  series,  leads  us  to  suspect  that  model  II  has 
a  drop  in  attenuation  which  is  too  large  between  250  and  700 
km  of  depth  (Fig.4). 


5. 


On  the  basis  of  the  results  shown  in  Fig.1,  we  infer  that  the 
above  observations  are  little  affected  by  the  finite  number 
of  radial  modes  used  in  this  study  since  the  phase  SS  appears 
to  be  converging  both  to  a  definitive  wavelet  form  and  ampli¬ 
tude. 

A  final  note  in  this  series,  which  will  describe  theoretical 
seismograms  generated  by  even  more  radial  modes,  will  give 
more  information  on  this  point. 

To  investigate  further  the  effect  produced  by  changes  in  the 
average  anelasticity  of  the  crust  and  upper  mantle,  we  gene¬ 
rated  time  series  (Fig. 3, trace  C)  corresponding  to  an  anela¬ 
sticity  distribution  derived  from  model  I,  but  with  the  pha¬ 
se  attenuation  B2  ,  reduced  to  0.000700  sec/lcm  from  the  top 
of  the  low-velocity  channel  down  to  the  400  Ion  discontinuity. 
Comparison  of  traces  C  and  A  shows  an  appreciable  decrease 
in  the  amplitude  ratio  S/Sa,  and  a  significant  increase  of 
the  energy  in  the  tail  of  the  Sa  wave  train. 

As  a  limiting  case  for  the  anelasticity  tests  we  have  also 
exhibited  the  time  series  (Fig.  3,  trace  D)  for  a  perfectly 
elastic  Earth.  In  this  case  ,for  reasons  reported  in  Part  II, 
the  surface  waves  completely  dominate  those  body  waves  that 
penetrate  deeply. 

The  phases  SS  show  a  less  drastic  decrease  in  amplitude  be¬ 
cause  they  are  relatively  strongly  affected  by  the  high  atte¬ 
nuation  in  the  upper  mantle. 

To  give  a  representation  of  the  effect  ,  on  the  time  series, 
produced  by  variations  in  source  depth,  we  present  theoreti¬ 
cal  seismoerams  (FI*?.  ^  - --  v  v  n  u\  ♦»<«* 
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res  as  traces  A,B,C,D  but  with  a  focal  depth  of  140  km. 

We  see  that  variation  of  the  source  depth  has  an  effect 
on  the  arrival  times  and  the  amplitudes  of  body  waves;  the 
feature  of  the  time  series, which  is  most  sensitive  to  the 
source  depth,  at  least  for  the  focal  and  structural  parame 
ters  we  used,  is  the  amplitude  of  the  SS  phase. 
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Fig. 


Fig.  2 


Fig.  3 


FIGURE  CAPTIONS 


Representation  of  the  effect  of  adding  more  higher  (radial) 
modes  in  the  synthesis  of  theoretical  sismograms  recorded 
at  the  surface  of  the  Earth.  Source  mechanism  is  a  dip 
slip  displacement  dislocation  on  a  vertical  fault  at  a  depth  T 
of  180  km.  Epicenter-station  separation  is  6788  km.  The 
model  we  have  used  for  the  intrinsic  anelasticity  of  the 
Earth  is  an  adaptation  (Fig. 4)  of  model  MM8  (Anderson,  Ben 
Menahem  &  Archambeau  ,  1965)  to  our  oceanic  structure. 

Torsional -wave  response,  at  the  surface  of  the  Earth,  to 
a  dip-slip  displacement  dislocation  on  a  vertical  fault, 
at  a  depth  of  180  km.  The  expected  body-wave  arrival  times 
are  indicated  for  s(  •  ),  SS  (A),  ScS  (▼),  sS  (■),  sscS 
(♦)  .  See  caption  of  Fig.  1  for  anelasticity  model.  ■. 

Comparison  of  theoretical  and  experimental  seismograms  ^ 

for  30-100  WWSSN  instruments,  6788  km  from  the  epicenter;  ~ 

short  traces  are  copies  of  the  long-period  experimental 

record  (SH)  measured  at  Honiara  (HNR)  for  the  event  occur-  - 

ring  atthefoot  of  the  Kamchatka  peninsula  at  14:30:30.3 

GCT  on  1964  December  26  (See  Kausel  et  al.,  1977  for  full 

details  on  this  event,  our  theoretical  treatment  of  the 

source  mechanism,  and  the  analysis  of  the  sa  portion  of  ^ 

the  time  series).  Traces  denoted  by  capital  letters  are  u 

theoretical  seismograms:  trace  A  represents  the  anelasti  »■ 

city  Model  I  (  Fig.  4  )  ;  trace  B,  anelasticity  Model  II  r~ 

(Fig.  4);  trace  C,  Ma-del  I  with  the  phase  attenuation  b2 

reduced  to  0.000700  sec/km  from  the  top  of  the  low  velo-  ! 

,  ! 
city  channel  down  to  the  "400-km"  discontinuity;  trace  D  L. 

does  not  Include  at-format-  1  nvi  f-  Vi  r»  t?  V*  *•»*>/*  *TVv»  s  I 


9 


E,  F,  G,  and  H  represent  the  same  models  respectively  as 
A,  B,  C  and  D  but  with  the  source  depth  reduced  from  160 
to  140  km. 

Pig.  4  S-wave  velocity  structure,  transverse-wave  phase  attenua¬ 
tion  B2,  and  intrinsic  attenuation  Q-1  ,  obtained  from  the 
adaptation,  to  our  velocity  structure,  of  two  anelasticity 
models  reported  in  the  literature:  (I)  world-wide  average 
model  (Anderson,  Ben  Menahom,  Archambeau,  1965),  (II)  0- 
ceanic  model  (Mitchell ,1 976 ) .  Model  I  is  given  by  solid 
lines.  Model  II  by  dotted  lines  . 


TABLE  1 

Amplitude-ratios  S/S  and  SS/S  measured,  peak  to  peak,  on  the 

cl  di 

traces  A,  B,  C,  E,  F  and  G  in  Fig.  3. 

1 80  km  1 40  km 


A 

B 

C 

E 

G 

S/Sa 

0.66 

0.62 

0.54 

0.59 

0.50 

0.43 

SS/Sa 

0.27 

0.39 

0.25 

0.48 

0.58 

0.41 
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Anderson 


Anderson 


Kausel  E., 


Mantovani 


Mitchell  B 


Nakanishi 


Schwab  F. 
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DEFINITIONS  OF  QUANTITIES  USED  IN  THIS  PAPER 


G 

u.'. 

U-'z 

fev 

Va, 

\ 

TU 

<T* 

•J2. 

CT2 

M, 

H 

K 

Q 

K 


Angular  deflection  of  seismometer  boom  in  radians 
Angular  deflection  of  galvanometer  mirror  in  radians 

Angular  frequency  of  seismometer  free  period,  =  ).rr  j  |t  ,  (  rad  i  ans /sec ) 

Angular  frequency  of  galvanometer  free  period,  =  H/  I  ^  ,  (radians/sec) 

Seismometer  damping  parameter,  =  Va,*  u_'i 

Galvanometer  damping  parameter,  =V\icu:i 

Seismometer  damping  constant  (1.0  =  critical  damping) 

Galvanometer  damping  constant  (1.0  =  critical  damping) 

Seismometer  coupling  constant 
Galvanometer  coupling  constant 

=  'T/'j’z  ,  coupling  constant  between  seismometer  and  galvanometer 
(1.0  =  d i rect  coup  1 ed ) 

Seismometer  boom  mass  (kilograms) 

Distance  between  center  of  gravity  of  mass  and  rotation  point 
of  seismometer  boom  (meters) 

Moment  of  inertia  of  seismometer  boom  (kg-rn~) 

Displacement  of  the  ground  (meters) 

=  Gu/H,  ,  magnitude  of  the  force  impressed  on  the  boom  by 
the  calibration  current  (Newtons) 


Calibration  current  (mi  I  I i amperes ) 

Q  Motor  constant  of  the  calibration  coil  (Newtons/Ampere) 

Xi’n  Height  of  the  impulse  response  as  measured  on  the  record  (milli- 
1  meters) 

V\( i ^  Heaviside  step  function  of  time 
HO  Dirac  delta  function  of  time 

. —  f  \  I  l  ( 

^  V\,  ,\\?  Pseudoparameter  forms  of  the  seismometer  and  galvanome- 
’  ter  free  periods  and  dampings 


The  time  at  which  the  calibration  current  was  switched  on 

An  apparent  amplitude  factor  that  is  determined  in  the  inversion 

A  constant  term  to  be  subtracted  from  the  digitized  data  to 
best  agreement  with  theoretical  impulse  responses 

A  linear  trend  to  be  subtracted  from  the  digitized  data  to 
give  best  agreement  with  theoretical  impulse  responses 

A  rotation  to  be  applied  to  the  impulse  response  to  correct 
for  incorrect  digitizing  of  the  impulse 


IMPULSE  RESPONSE  INVERSION 


Quantitative  time  series  analysis  using  the  World-Wide  seismic 
network  records  (WWSSN)  requires  that  corrections  be  made  for  the 
response  of  the  seismograph.  In  the  usual  method  of  calculation  of 
instrumental  corrections  the  seismograph  is  considered  as  a  black  box. 
From  Fourier  analysis  of  the  response  to  a  known  transient  impulse 
(usually  a  step  function  in  acceleration  applied  to  the  seismometer 
mass)  corrections  can  be  obtained.  In  general,  at  frequencies  in  the 
neighborhood  of  the  peak  instrumental  response  this  is  an  accurate 
ethod  of  correction.  However,  periods  shorter  than  the  seismometer 
free  period  are  not  strongly  excited  in  the  impulse  response  and  sig¬ 
nificant  errors  in  the  computation  of  ground  motion  may  occur. 

Studies  involving  short  period  crustal  and  higher  mode  surface 
waves  have  obliged  us  to  seek  an  alternative  method  of  determining  these 
corrections,  in  order  to  reduce  errors  caused  by  instrumental  uncer¬ 
tainties.  The  parameters  of  the  seismometer  system  have  been  determined 
directly  using  a  linear  inversion  method.  Mitchell  and  Landisman's 
(1969)  analytical  expressions  can  then  be  used  to  derive  the  phase  and 
amplitude  corrections  at  all  periods.  This  method  uses  a  stabilized 
I eas t- squa res  technique  which  allows  more  parameters  to  be  used  while 
maintaining  convergence  in  the  inversion  process.  Arbitrary  values  of 
the  coupling  constant  and  the  seismometer  period  can  be  used  in  this 
form  of  the  calculation.  Also,  the  scale  factor  determined  in  the  in¬ 
version  can  be  used  to  obtain  an  expression  for  seismograph  magnifica¬ 
tion  that  is  independent  of  the  rated  value. 
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The  procedure  for  inversion  ( J  acksori ,  I  9  72  )  involves  solving  a  sys 
tern  of  linear  equations  for  the  corrections  Ax'  "to  the  model  para- 
jters.  In  our  case  these  parameters  are  the  seismograph  instrumental 
constants.  The  equations  are  set  equal  to  the  data  residuals  j 

which  are  the  residuals  between  the  data  and  the  predictions  from 
the  model.  The  linear  equations  take  the  form 


-  [~\i^  /X  1*  £.1 

where  the  matrix  ‘  .  tfle  partial  derivatives  of  the 

data  residuals  with  respect  to  the  model  parameters.  The  1 eas t- squa res 


solution  results  from  the  s t ra i ght forward  minimization  of  the  error 
siduals  (Mitchell  and  Landisman,  1969)-  They  noted  that  the  in¬ 
clusion  of  too  many  parameters  tended  to  cause  instabilities  in  the 


inversion.  Thus  certain  parameters,  such  as  the  scale  factor  and 
impulse  'on'  time,  had  to  be  empirically  adjusted  in  the  inversion. 

In  our  experience  with  digitizing  actual  impulse  responses  there  is 
error  introduced  due  to  baseline  choice,  and  so  a  constant,  linear 
2nd,  and  rotation  parameters  were  introduced.  Although  all  the 
partial  derivatives  have  linearly  independent  components,  in  the 
presence  of  noise  not  all  of  them  can  be  resolved,  and  an  unstable 
inversion  resulted  using  I eas t- squares  method  with  the  nine  parameters 
listed  in  figure  I.  The  inversion  was  stabilized  by  decreasing  the 
degrees  of  freedom  in  the  inversion  (Wiggins,  1972).  The  matrix 

of  partial  derivatives  is  first  normalized  with  respect  to  the  data 
standard  deviations  and  _a  p r i o r i  uncer ta i n t i es  assigned  to  the  model 

ameters.  The  matrix,  now  in  normal  form  (Jackson,  1972),  is  factored 
according  to  the  Lanczos  f actor i zat i on  (Lanczos  1961)  into  eigenvalue 
and  eigenvector  matrices.  In  this  form  eigenvalues  less  than  1.0 
correspond  in  some  sense  to  parameter  combinations  that  cannot  be 


-3- 


resolved  from  the  data,  and  so  th-.se  eigenvalues  are  set  to  zero, 
removing  these  degrees  of  greedom  from  the  inversion.  In  general 
the  number  of  independent  degrees  of  freedom  in  the  inversion 
when  this  was  done  was  six  or  seven.  All  the  parameters  necessary 
to  obtain  phase  and  magnification  information  were  obtainable, 
but  careful  digitization  was  necessary  to  insure  good  results. 

COMPUTATIONS 

The  actual  computation  of  the  partial  derivatives  was  done 
using  an  analytical  formula  for  the  inpulse  response  as  a  func¬ 
tion  of  the  parameters.  Mitchell  and  Land i srnan ( 1 969 )  used  the 
inverse  Fourier  transform  of  the  frequency  domain  representation 
of  the  impulse  response  to  calculate  the  partial  derivatives. 

While  this  has  certain  advantages,  I  felt  that  the  number  of  Fou¬ 
rier  transforms  needed  per  inversion  would  be  a  significant  cost 
factor  even  with  FFT,  and  so  a  method  was  followed  as  illustrated 
in  figure  2.  Starting  with  the  representation  of  the  seismometer 
and  galvanometer  as  coupled  damped  harmonic  oscillators,  the 
equation  of  motion  of  the  rotation  of  the  galvanometer  is  solved 
by  inverse  Laplace  transform.  The  magnitude  of  the  impulse  given  to 
the  mass,  ,  is  just  •  Jarosch  and  Curtis  (1C' ’’3)  noted 

that  the  denominator  cannot  be  factored  unless  the  term  containing 
(j'i  is  set  to  zero.  When  this  is  done  the  pseudo-parameter 
equation  can  be  solved  uniquely  to  obtain  new  values  of  ,  ,1 /_  , 

and  we  call  the  new  values  of  I,  (lj  ,\n,  *  the  1  pseudoparameters  1  . 

As  can  be  seen  from  the  figure,  the  effect  of  this  conversion 
is  small  for  <XZ  small,  but  may  be  significant  for  larger  r)  *  . 
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Furthermore,  the  pseudoparameters  become  complex  for  larger 

We  can  obtain  finally  an  analytic  expression  for  the  impulse  response: 


where 

This  will  be  in  general  comp  1  ex- va 1 ued  for  complex  parameters; 
the  real  part  corresponds  to  the  impulse  response  function. 

If  the  seismometer  or  galvanometer  becomes  critically  damped, 
we  must  remember  that  K\rO  or  Ki  -  0  ,  and  so 

S-^7~-  "  i  ^  •  CC‘S  1  }  COS  Kit  - 


For  the  case  when  the  seismometer  or  galvanometer  is  overdamped, 
or  ^  becomes  complex,  and  so  the  appropriate  sin  and 
cos  are  replaced  by  sinh  and  cosh.  In  figure  3  are  some  impulse 
responses  calculated  using  the  above  formula.  The  reference  impulse 
response  corresponds  in  some  way  to  a  standard  WWSSN  instrument. 

Some  notion  of  what  the  partial  derivatives  look  like  can  be 
gained  from  this  figure.  The  effect  of  changing  the  seismometer 
free  period  can  be  seen  to  be  very  similar  to  the  effect  of  chang¬ 
ing  the  gain,  or  amplitude  factor  of  the  impulse  response.  This 
effect  can  be  seen  in  the  Lookout  Mountain  calibration  tests,  in 
which  two  different  impulse  responses  from  the  instrument  were 
inverted  independently.  The  seismometer  free  period  values  show 
somewhat  larger  discrepancies  than  the  other  parameters.  The  theo¬ 
retical  values  for  the  parameters  compare  very  well  with  those 
determined  experimental  I y  using  the  usual  tests  for  damping  and 
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free  period.  In  the  case  of  the  galvanometer  constants,  there  is 
a  very  good  argument  for  determining  the  constants  through  inversion 
rather  than  by  other  methods,  since  this  seems  to  be  at  least 
as  accurate  and  does  not  disturb  the  instrument  in  any  way. 

The  seismometer  constants,  especially  the  free  period,  are  probab¬ 
ly  better  determined  by  other  experimental  methods,  and  constrained 
in  the  inversion.  Once  the  instrumental  constants  are  determined, 
the  instrumental  effects  on  amplitude  and  phase  can  be  determined 
theoretically.  Following  Mitchell  and  Landisman  (.1969)  the  phase 

can  be  w  r i t  ten  as : 

I//"-  )  <  -V  h,',j  r  4l  (l- tr2')]  to'  c 

t/v  i  CU_)  —  v  ,  ,  .  ,  £  ,  , 

■  A u>'*  —  A(  )  oj 

The  convention  here  is  that  c^/i^u)  is  positive,  which  corres¬ 
ponds  to  positive  instrumental  group  delay  time.  The  magnification 
of  the  instrument  can  similarly  be  expressed  in  terms  of  the  para¬ 
meters  of  the  instrument: 

I'lMwV’/G  if) 


I  V  >  >  i 

U,  I  u.  ; 


is  the  inversion  gain  constant,  ip  is  the  calibration  coil 
current,  and  Q  is  the  motor  constant  of  the  coil.  For  WWSSN 
stations  there  is  an  empirical  formula  for  the  gain  at  the  seis¬ 
mometer  free  period: 

M~  X(.p  for  T,  =15  seconds; 

Ca  Ip 


H- 


_  P* 1  ^ 


Gif 


for  Ti  =30  seconds 


■^■^is  the  height  of  the  impulse  response  in  mfl 1 imeters .  This 


formula  was  checked  against  the  analytical  formula  for  gain  and 
agreement  to  a  few  percent  or  better  was  found.  However,  the 


i 
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empirica!  formula  assumes  all  the  constants  are  at  the  rated  values 
and  so  a  badly  tuned  instrument  would  give  erroneous  gains. 

A  long  period  WWSSN  station,  NDI-LPZ,  was  used  in  inversion  to 
get  instrumental  parameters.  The  amplitudes  in  the  impulse  response 
must  be  multiplied  by  UJ'>  to  get  the  ground  magnification.  This 
means  that  the  gains  at  the  long  period  end  are  well  determined, 
where  the  original  impulse  had  a  strong  signal.  However,  below 
30  seconds  period  the  amplitudes  become  relatively  worthless  and 
it  is  advisable  to  use  the  theoretical  values.  The  same  consider¬ 
ations  hold  for  the  phases.  An  additional  source  of  error  In  the 
phase  of  the  impulse  is  introduced  because  the  'on'  time  of  the 
impulse  is  not  known  a  priori ,  but  must  be  estimated  from  the  re¬ 
cord.  An  error  of  one  second  or  more  is  common  here.  This  leads 
to  a  systematic  bias  in  the  phases  determined  as  can  be  seen  in 
figure  5.  where  the  empirically  determined  'on'  time  is  wrong 
by  one  second.  No  bias  is  evident  in  the  phases  derived  using 
the  'on'  time  from  inversion. 

A  note  concerning  baseline  determination  is  relevant  here. 
Although  it  was  not  possible  in  the  inversion  to  use  the  rotation 
as  an  independent  parameter,  the  rotation  is  quite  important,  es¬ 
pecially  for  very  large,  steep  impulses,  and  a  wrong  baseline  will 
distort  the  constants  obtained  and  resultant  phases  quite  badly. 
However,  drawing  the  correct  baseline  is  not  juu  a  matter  of  con¬ 
necting  the  ends  of  the  trace  and  making  sure  the  impulse  does  not 
occur  in  a  disturbed  part  of  the  record.  Figure  7  shows  an  impulse 
on  our  instrument.  The  minute  marks  line  up,  which  makes  it  easy 
to  notice  the  offset  of  the  minute  mark  near  the  top  of  the  pulse. 
This  is  probably  due  to  a  lens  misalignment.  This  effect  may  be 
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common  among  WWSSN  instruments,  as  noted  by  James  and  Linde  (1971)- 
What  I  did  for  this  was  to  rotate  the  impulse  on  the  digitizer  un¬ 
til  the  minute  marks  on  the  pulse  fell  at  equal  intervals  on  the 
digitizer  scale.  This  meant  that  the  digitizer  scale  was  then 
aligned  parallel  to  the  direction  of  motion  of  the  galvanometer 
light  beam.  Digitizing  the  impulse  in  this  manner  gave  satisfactory 
results  even  for  large  impulse  responses. 

Some  conclusions  I  have  reached  in  this  work  are  that  the  appli¬ 
cation  of  a  generalized  inversion  scheme  to  instrument  impulse 
responses  allows  in  inclusion  of  more  parameters  without  sacrifi¬ 
cing  the  stability  and  validity  of  the  inversion.  In  particular, 
the  use  of  the  amplitude  factor  allows  an  independent  measure  of  the 
magnification  to  be  macle.  The  inclusion  of  the  'on'  time  means 
that  this  source  of  phase  error  can  be  practically  eliminated. 

The  constant  and  linear  trend  terms  are  important  because  of  digi¬ 
tizing  cons i derat  ions  in  which  a  baseline  must  be  assumed  and  later 
corrected  for.  Careful  digitizing  is  necessary  to  assure  that  the 
correct  baseline  is  used. 

Noise  on  the  record  during  the  impulse  response  did  not  seem 
to  be  a  significant  facto  i  r»  the  inversion.  In  general  I  chose 
only  impulses  that  occured  during  clear,  quiet  times  and  without 
visible  noise  except  for  microseisms.  Examining  the  residuals  that 
remain  after  inversion,  I  found  that  much  of  the  error  seemed  to  be 
associated  with  digitizing  errors  on  the  steep  parts  of  the  impulse. 
This  is  unavoidable  because  of  the  steepness  of  the  pulse  and  the 
uncertainty  in  the  digitizer  scale.  The  largest  magnitude  of  these 
errors  appeared  to  be  of  the  order  of  1%  of  the  impulse  height 
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for  good  impulses  that  had  successful  inversions.  These 
errors  did  not  appear  to  cause  much  error  in  the  parameters 
determined  as  shown  by  the  good  agreement  between  indepen¬ 
dent  impulse  inversions  from  the  same  instrument. 

ERROR  ANALYSIS  AND  CONCLUSIONS 

Errors  in  the  determination  of  parameter  values  are  of 
several  different  types-  One  type  of  error  arises  because 
of  the  finite  line  thickness  and  other  digitizing  errors,  and 
because  of  microseisms  and  other  types  of  noise  present  on 
the  record.  These  effects  can  all  be  included  as  noise  in 
the  impulse  response  signal.  Figure  8  shows  the  results  of 
inversion  of  theoretical  and  actual  impulse  responses. 

Random  noise  (white  noise)  was  added  to  one  theoretical 
impulse  and  inversion  was  done.  Convergence  of  error  to  the 
noise  level  was  quite  rapid,  and  after  that  essentially  no 
change  took  place.  Similar  results  were  obtained  for  two 
actual  impulses,  LMO-2  and  NDI-LPZ.  The  digitizing  noise 
level  was  obtained  by  an  estimation  of  the  random  error  to 
be  about  .  1 4  mm  on  the  record  on  the  average.  However,  in 
the  steep  parts  of  the  impulse  the  errors  became  greater,  and 
residuals  of  amplitude  equal  to  about  1%  of  the  impulse 
response  amplitude  (about  I  mm)  were  occasionally  seen. 

Since  the  steep  regions  of  the  impulse  response  are  most 
critical  in  determi nat i on  of  the  instrumental  parameters, 
this  upper  bound  was  taken  as  the  representative  noise  level 
in  the  analysis  of  propagation  of  errors.  The  effect  of 
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random  errors  in  the  data  on  the  parameters  was  estimated 
by  computing  the  covariances  of  the  model  parameters  by  the 
method  of  Jackson  ( 1 9 72 )C 

oav  (  t:  l_,  V\v\ 

The  matrix  \A  is  r  \J  /\  \J  ,  the  general  inverse 
of  the  matrix  of  partial  derivatives.  The  results  obtained 
using  I  mm  as  the  random  error  are  shown  in  table  1.  The  in¬ 
strumental  parameters  for  LMO- 1  and  LMO-2  appear  to  agree 
within  the  uncertainties  given  for  the  parameters,  although 
the  errors  in  the  data  may  not  actually  be  random  as  assumed. 

The  errors  on  the  seismometer  and  galvanometer  free  periods  are 
about  equal  for  the  severely  underdamped  case,  LMO- 1  and  LMO-2, 
but  change  for  the  overdamped  case,  LMO-NEW.  The  impulse  response 
was  activated  by  an  automatic  relay  at  exactly  12-hour  intervals. 
There  appears  to  be  some  difference  between  the  'on'  times  of  the 
impulses  LMO- 1  and  LMO-2  on  the  seismogram,  as  judged  by  measure¬ 
ments  relative  to  the  minute  marks.  This  may  account  for  the 
difference  between  the  "T0  times  of  LMO- 1  and  LMO-2  impulse 
responses.  Q  ,  the  rotation,  was  constrained  to  equal 
zero  throughout  the  inversions. 

Another  type  of  error  is  not  random,  and  results  because  the 
seismometer  free  period  and  the  amplitude  scale  parameter  p\  can¬ 
not  be  resolved  from  the  data.  Because  of  this  trade-off  effect, 
errors  due  to  lack  of  resolution  will  have  little  effect  on  the 
calculated  magnifications  or  phases.  It  could  be  corrected  by 
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constraining  in  the  inversion,  since  this  parameter  can 

be  accurately  measured.  Table  2  shows  the  result  of  constraining 
Xt  by  adding  in  the  parameter  eigenvector  to  the  solution. 

X  and  Pi  have  both  been  made  consistent  by  this  procedure. 

Some  conclusions  that  can  be  drawn  from  this  work  are  that  ap¬ 
plication  of  generalized  linear  inverion  to  seismograph  impulse 
responses  to  derive  the  instrument  constants  yields  useful  and 
consistent  results,  and  makes  possible  more  accurate  estimates 
of  gain  and  phase  shift  at  short  periods.  An  independent  means 
of  deriving  the  gain  as  a  function  of  period  is  possible  using 
the  inversion  amplitude  parameter.  Errors  in  phase  due  to  un¬ 
certainty  in  impulse  'on'  time  are  reduced.  There  is  a  trade¬ 
off  of  ,  the  seismometer  free  period  and  f\  ,  the  gain 

constant,  but  these  can  be  resolved  by  constraining  X  .  which 
is  in  general  well-known.  Noise  caused  by  digitizing  errors  or 
microseisms  does  not  appear  to  be  a  problem  in  this  work.  The 
residuals  were  usually  of  amplitude  less  than  1%  of  the  original 
impulse  amplitude.  It  was  necessary  to  take  great  care  in  choos¬ 
ing  the  correct  baseline  for  digitizing.  Application  of  the 
inversion  technique  to  actual  impulse  responses  yielded  consis¬ 
tent  results  and  enabled  amplitude  and  phase  shifts  to  be  cal¬ 
culated  at  all  periods. 


BIBLIOGRAPHY 


Burdick, L.  and  G.  Mellman,  The  response  of  Che  WWS5N  short 
period  seismometer  (in  press). 

Espinosa,  A.F.,  G.H.  Sutton  and  H.J.  Miller,  S.J.,  A  tran¬ 
sient  technique  for  seismograph  calibration  -  manual  and  stan¬ 
dard  set  of  theoretical  transient  responses,  VESIAC  special  re¬ 
port  No.  4410- 106X. 

Hagiwara,  T-  (1958),  A  note  on  the  theory  of  the  electro¬ 
magnetic  seismograph,  Bull.  Earthquake  Res.  Inst.  _3§,  139-164. 

Jackson,  D.D.  (1972),  Interpretation  of  inaccurate,  insuf- 
fient  and  inconsistent  data,  Geophys.  J.R.  astr.  Soc.  _23,  97-110. 

James,  D.E.,  and  A.T.  Linde  (1971),  A  source  of  major  error 
in  the  digital  analysis  of  World-Wide  Standard  Station  seismo¬ 
grams,  Bull.  Seis.  Soc.  Am.  _6]_,  723-  728. 

Jarosch,  H.  and  A. R.  curtis  (1973),  A  note  on  the  calibra¬ 
tion  of  the  electromagnetic  seismograph,  Bull,  Seis.  Soc.  Amer. 
63,  1145-1155- 

Lanczos,  C.  (1961),  Linear  differential  operators,  D.Van 
Nostrand  Co.,  London. 

Mitchell,  B.J.  and  M.  Landisman  (1969),  Electromagnetic 
seismograph  constants  by  I eas t- squares  inversion,  Bull.  Seis. 

Soc.  Am.  59.  1335- 1348. 

Operation  and  maintenance  manual  World-Wide  Seismograph 
System,  model  10700,  The  Geotechnical  Corporation,  3401  Shiloh 
Road,  Garland,  Texas. 

Wiggins,  R.A.  (1972),  The  general  linear  inverse  problem: 
implication  of  surface  waves  and  free  oscillations  for  Earth 
structure,  Rev.  Geophys.  Space  Phys.  j_0,  251-285- 


APPENDIX  A 


ANALYTICAL  METHOD  FOR  OBTAINING  PSEUD CPA RAMETERS 


From  figure  2  we  have  Che  equality  involving  original  and 
pseudoparameters  that  must  be  satisfied: 

Following  Jarosch  and  Curtis  (1971)  we  equate  powers  of  S 


(1)  €/+?-'  = 

(2)  ‘t-E.'tV  -  -  Mf.f.r2 

o)  u,\,a*;  ■+u>;i‘c.'  - 

(4)  u.\,2w':z-  u2u'* 

To  simplify,  let 

(5)  W?  -  cl 

(6)  u':2ui.;se 

(7) 

(8)  c\'t>'4iT3'  '  c 

(9)  L-' J E 4  -v-  Vwj.'**!  - 

We  can  t:  derive  the  cha rac ter  I s c i c  equation  in  CX 

(10)  -'CO?  ^  •"  c?  -  1 1  (’t  —  cl  -c4l  )  cxJ  CM  W4  -  |''V<X  — ( 

Rather  than  solve  this  equation  numerically  as  Jarosch  and  Cur¬ 
tis  did,  I  needed  the  complex  solution'  tc  the  equation.  The 
equation  can  be  factored  as 

(11)  (  a*-*  (o*  +^<s  -'l 


Si nee  there  is  symmetry  between  x  ,  l 


We  must  solve  for  x,  j  ,  £ 
they  may  be  interchanged; 


(12)  X,  Vj,-L  ~  (\  *8 


+  (t-VQK^ 

2.  2 


- ')  -  (h  -eo  PT 

"-JL  2 


where  * 

( :  j)  f-\  -  t  S)  S 

(14)  8  -  (-b/;  -  pTPt//27'  P-3 


and 

(15)  ft  -.8 

(16)  Vj  r  2n'J4'?  ~  |3r/S  -  °< 

(17)  c*  -  4  tp  ■- *W2- 

(18)  /}  ~  *  p 

The  solution  for  10**"  is  _ ...  , 

i0<i;  t'*~’ZrvN{,z 

(19)  5  1,-xj  i  VP 

~  [-1  * 

I  always  used  the  positive  sign  in  the  root.  If  one  root  was  t0>' 
the  other  two  corresponded  tot^  and  u>..  It  was  always  possible 
to  distinguish  which  was  which. 


FIGURE  CAPTIONS 


Figure  I:  Schematic  mode)  of  the  impulse  response  inversion 
program.  The  periods  and  dampings  are  in  pseudopa remoter  form 
as  used  in  the  inversion.  Up  to  five  iterations  were  used  in 
convergence  to  the  stabilized  1 eas t- squa res  solution. 

Figure  2:  Equations  for  the  seismometer  and  galvanometer  motions, 
following  Jarosch  and  Curtis  (1973)-  The  effect  of  the  pseudo¬ 
parameter  conversion  is  small  for  small,  but  significant 

for  larger 

Figure  3:  Effect  of  various  parameters  on  the  impulse  response 
function.  The  solid  line  is  the  for  all  figures.  Changing  the 
seismometer  free  period  has  a  large  effect  on  the  impulse  am- 
pl i tude. 

Figure  4:  Results  of  inversion  on  our  long  period  vertical  in¬ 
strument.  The  parameters  were  measured  in  situ  using  overshoot 
ratio  tests  and  timing  the  free  periods  Fy  stopwatch. 

Figure  5:  Phase  determination  of  NDI-LPZ  impulse  response,  using 
direct  Fourier  transform  method  and  the  inversion  method.  The 
Fourier  method  is  degraded  at  periods  below  30  seconds.  The 
'on'  time  had  been  picked  by  eye  at  4  seconds  after  the  minute 
mark;  the  inversion  showed  it  actually  occured  at  3  seconds  after 
the  minute.  The  large  black  dots  correspond  to  the  best  esti¬ 
mate  of  the  phase. 

Figure  6;  Amplitude  determination  of  NDI-LPZ  using  the  Fourier 
method  and  inversion  method.  The  curve  on  the  left  corresponds  to 
impulse  an'.o  1  i  tudes .  These  become  very  small  below  30  seconds, 
which  accounts  for  the  degradation  of  the  phases  and  amplitudes. 

To  obtain  ground  displacement  magnification,  this  curve  is  mul¬ 
ti  p  I  i  ed  by 

Figure  7:  Impulse  response  on  the  Lookout  Mountain  (LMO)  LPZ 
record.  There  is  a  slight  misalignment  of  the  light  beam  rela¬ 
tive  to  the  drum  as  shown  by  the  offset  of  the  minute  mark  on 
the  impulse.  This  minute  mark  was  used  to  help  align  the  impulse 
for  correct  digitization. 

Figure  8:  Convergence  of  the  linear  inversion  method  for  theo¬ 
retical  and  actual  impulse  responses.  For  all  successful  solu¬ 
tions  the  method  converged  within  5  iterations. 
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COMPUTATION  OF  COMPLETE  THEORETICAL 
SEISMOGRAMS  FOR  TORSIONAL  WAVES 

By  A.  H.  Liao,  F.  Schwab,  and  E.  Mantovani 
Introduction 

In  Che  study  of  earth  structure,  and  earthquake  source  mechanism,  the  handiest 
and  probably  most  important  information  is  that  provided  by  the  seismic  waves 
recorded  on  seismograms.  Therefore,  the  direct  comparison  between  theoretical 
and  experimental  seismograms  is  a  very  appealing  subject  for  most  seismologists. 

The  computation  of  theoretical  seismograms  has  been  studied  by  various  authors. 
Some  tried  to  obtain  long-period  seismograms  by  combining  a  few  normal  nodes; 
some  tried  to  obtain  a  few  body-wave  phases.  The  results  of  the  investigations 
that  led  to  the  present  letter  are  reported  by  Schwab  (1970);  Schwab  and  Knopoff 
(1970,  1971,  1972.  1973);  Kausel  and  Schwab  (1973);  Schwab  and  Kausel  (1976); 
Knopoff  et  at.  (1973);  Knopoff  er  al.  (1974);  Schwab  et  al.  (1974);  Nakanishi  et  al. 
(1976);  Kausel  et  al.  (1977);  Nakanishi  et  al.  (1976);  Mantovani  et  al.  (1976); 
Mantovani  (1977a);  and  Mantovani  k  in  preparation).  The  results  of  our  most 
recent  work  now  allow  us  to  generate  complete  theoretical  seismogTams  for  torsional- 
wave  motion.  By  saying  complete,  we  mean  that  all  modes  that  exist  for  periods 
above  10  sec  are  included  in  our  seismograms,  i.e.,  that  all  amplitude  and  phase 
information  down  to  a  period  of  10  sec  is  included.  This  means  that  the  body-wave 
phases,  as  well  as  surface-wave  arrivals,  are  obtained  with  this  method. 

The  computational  technique  is  somewhat  involved;  however,  anyone  who  is 
interested  in  the  details  can  obtain  program  copies  from  the  authors.  Briefly,  we 
compute  all  of  the  required  frequency-domain  information  such  as  phase  and  group 
velocity  dispersion,  attenuation,  amplitude  excitation,  apparent  initial  phase,  depth 
of  penetration  (and  if  desired,  partial  derivatives  with  respect  to  structural  param¬ 
eters!  from  the  longest  existing  period  for  each  mode,  down  to  10  sec.  This  requires 
between  90  and  100  radial  modes,  and  is  accomplished  in  a  single,  relatively  short 
computer  run.  This  computation  corresponds  to  about  50,000  free  oscillating  modes. 
Theoretical  seismograms  in  the  time  domain  are  then  obtained  by  inverse  Fourier 
transformation  as  described  by  Kausel  and  Schwab  (1973)  and  Schwab  and  Kausel 
(1976).  The  cost  of  computing  the  volume  of  frequency-domain  information,  which 
is  required  to  construct  the  synthetic  seismograms  down  to  a  period  as  short  as  10 
sec,  has  been  prohibitively  high  up  to  the  present  time:  but,  with  our  current 
optimization  we  have  succeeded  ui  reducing  the  computation  time  to  about  6  min 
on  our  IBM  360/91  computer  (an  expenditure  of  about  $50)  for  a  given  earth 
structure.  The  corresponding  time  on  a  CDC  7600  computer  is  about  2  min.  [For 
corresponding  times  on  other  computers  see  Schwab  et  al.  (1977)  and  Porter  et  al. 
(in  preparation).]  These  time  estimates  should  be  considered  only  as  upper  bounds; 
relatively  simple  improvements — which  a  little  practical  experience  will  be  needed 
to  justify — are  expected  to  cut  these  times  at  least  in  half. 

The  purposes  of  the  present  letter  are  (1)  to  report  that  the  capability  now  exists 
for  computing  realistic,  torsional-wave  seismograms  that  contain  all  of  the  seismic 
energy  and  arrivals  down  to  a  period  of  10  sec;  (2)  to  announce  the  availability,  for 
general  distribution,  of  the  associated  program;  and  (3)  to. exhibit  the  results  of  our 
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initial  tests  of  direct  comparison  of  these  complete  theoretical  seismograms  with 
an  experimental  record  from  a  standard,  long-penod  instrument  of  the  WWSSN. 


Examples 

For  our  comparison  of  theoretical  and  experimental  seismogTams,  we  have  chosen 
an  event  with  its  epicenter  near  the  foot  of  the  Kamchatka  peninsula  [see  Kausel 
et  al.  (1977)  for  a  discussion  of  this  event).  Analysis  of  this  earthquake  indicates 
that  the  source  is  approximated  by  purely  dip-slip  motion  on  a  vertical  fault  plane. 
We  used  the  station  at  HNR,  which  provides  a  N-S  path,  so  that  the  recording  on 
the  E-W  instrument  can  be  considered  to  contain  only  the  azimuthal  component 
of  displacement.  The  focal  depth  is  in  the  intermediate  range.  /The  initial  earth 
model  that  we  used  in  our  calculations  was  a  somewhat  smoothed  version  of  the 
average  oceanic  structure  described  in  Kausel  et  al.  (1977),  which  is  essentially  the 
CIT-ll  structure  (Figure  10,  Toksoz  and  Anderson,  1966).  About  200  layers  were 
used  to  model  the  vertical  heterogeneity  of  the  crust-mantle  system;  an  adaptation 
of  model  MM8  (Anderson  et  al.  1965)  was  used  to  approximate  the  dependence  of 
the  intrinsic  attenuation  on  depth. 

Tk  prepare  the  experimental  and  theoretical  seismograms  for  comparison,  the 
usual  ramp  function  was  removed  from  the  experimental  time  series  (James  and 
Linde,  1971).  The  result  of  this  is  the  upper  trace  in  Figure  1.  Since  the  theoretical 
seismograms  do  not  contain  periods  below  10  sec,  we  ran  a  low-pass  frequency 
(boxcar)  filter,  with  cutoff  at  10  sec,  over  the  experimental  data  (second  trace  in 
Figure  1).  A  Gaussian  roll-off  filter  with  peak  at  15  sec  and  a  90  per  cent  decay  at 
10  sec,  was  then  applied  to  suppress  the  time-domain  oscillations  introduced  by 
the  abrupt  cutoff  of  the  boxcar  filter.  The  result  of  these  operations  on  the 
experimental  seismogram  is  shown  in  the  third  trace  of  Figure  1  The  same  Gaussian 
roll-off  filter  is  applied  to  the  theoretical  seismograms,  so  that  both  the  theoretical 
and  experimental  time  series  are  processed  in  the  same  way.  Thus,  we  can  compare 
these  time  series  directly. 

In  our  sample  computations,  we  have  varied  the  souice  parameters  over  the 
possible  range  of  variation  from  our  point-source  solution  The  first  figure  shows 
the  dip  angle  dependence.  With  the  theoretical  seismograms  available  for  compari¬ 
son,  we  can  easily  identify  several  of  the  arrivals  on  the  experimental  record.  The 
most  striking  of  these  are  the  body-wave  phases  S,  sS.  ScS,  sScS.  and  the  later, 
surface-wave  arrivals.  (The  various  arrivals  are  identified  in  Figure  5.)  For  this 
particular  event,  the  surface  waves  are  composed  of  multimode,  oceanic  So  (Kausel 
et  al.,  1977).  In  Figure  1  we  have  varied  the  dip  angle  from  that  of  a  vertical  (90°) 
fault  plane,  to  70°:  the  results  show  that  as  dip  angle  decreases,  the  normalized  S 
and  ScS  phases  remain  unchanged,  but  sS  and  sScS  dimmish.  Thus,  the  dip  angle 
appears  to  have  its  main  effect  on  body  waves  when  they  leave  the  focus  above 
the  equator  of  the  focal  sphere.  Both  the  long-  and  short-period  energy  increase,  m 
the  surface-wave  portion  of  the  seismogram,  when  the  dip  angle  decreases  from 
vertical.  Figure  2  illustrates  the  dependence  of  the  seismogram  on  source  depth. 
Here  we  have  varied  this  depth  from  220  to  140  km:  the  consequent  separation  of 
bodv-wave  arrivals,  and  development  of  the  Sa  wave  train  in  the  lower  set  of 
theoretical  traces,  clearly  exhibit  the  power  of  complete  theoretical  seismograms 
as  interpretive  tools.  As  Figure  3  illustrates,  the  theoretical  seismograms  from  this 
class  of  sources  are  less  sensitive  to  changes'in  slip  angle  than  to  changes  in  fault 
dip  or  source  depth. 

The  azimuthal-component  seismogram  is  composed,  in  principle,  of  two  contri- 
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butions:  that  from  the  torsional-mode  ejtcitation.  and  that  from,  the  spheroidal  ,  ; 
modes.  The  general  depressions  for  the  polar  \§)  and  azimuthal  ($)  components  of  * 
displacement  are  given  below  for  the  torsional  (D  and  spheroidal  iS)  excitations 
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FlO.  1.  Dependence  of  computed  *ei>moeram*  on  dip  ancle  it  ’He  -ow  -  For  r.„,,raH,  •*n»*r» 

mental  trace*  ettend  further  in  th*-  left  and  nch'.  and  tb*  dominant  le.iMward'  M>e  of  the  •ir«f 
amval— S—  u  upward,  on  the  *h»rter  theoretical  seiamocTam*.  •■luuwa/d  i*  down  Ka»  h  trn<  e  *rn.>fm, 
ued  relative  to  it*  maximum  displacement  from  rem  From  th*  »on  down  tli-  rip*'tn<.<niat  ir*.  r,  ,h-  „ 
the  K-W  lone -period  seismogram  at  HNH  with  the  ramp  function  removed,  then  with  a  low  pass 
frequency  filter  applied  icutolf  at  0  I  Hn  alter  ramp  removal,  and  the  l.mer  4^  hKv«.  a  <,,ius«,.,n 
roll  off  filter  applied  In  Figure*  1  to  1  we  show  the  effect.*  of  va/vtng  aource  parameter*  from  an  initial 
fm-  •  m  ancle  n*  *h*  fault  plan,-  h-low  »h-  hnm-.nt.al.  4  nf  <*>"  •nmuthml  ancle  between 

•tnke  line  and  epwrenler-to-ai  it  ion  line.  *.  of  40*.  wMtrt*  depth,  h.  i*  uutinllv  set  at  lr*0  km  ihau«el  er 
of  I9T7.  give  the  rra*on  for  thi*  initial  hmre  nf  source  dept  hi.  dip  angle  of  the  hanging  well  ret  mve  tn 
the  foot  wall.  \.  of  9»»°.  and.  a  step-hinclion  dieplacement  dislocation  is  gummed  at  the  focus.  *hert-  * 
point  aource  i*  avmmed.  v  well  a*  a  component  nf  *tres«.  normal  tn  the  fault  plane,  which  t*  contimnHi* 
dunng  the  dialocatmn  fSee  Figure  h*n*e|  et  ul  •  IS  r  ‘  i  |..r  fault  pi  me  c-o-*«.irv  j  In  »he  ,»  t%k?1»^ 

the  awHinted  dip  anfle  l*  civen  at  the  left  of  each  theoretical  tru  e  the  upper  .he..-  i,,,,,. 

sene*  represents  actual  ground  displacement  and  the  lower  *et  incindr*  in.'mmrnui  response  (the  u» 
("to  100  instrument  desrnbed  by  Knu|Mfl  tt  al.  1 1973) |  and  Cau.vu«n  roil  utf  filter 
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with  m  S  2.  Since  (u4)r  and  (u«).s  are  of  the  same  order  of  magnitude,  so  must  Vfr 
and  Vs.  Thus,  in  terms  of  ail  order-of-magnitude  estimate,  the  relative  sizes  of  the 
two  contributors  to  azimuthal-component  seismograms.  (uyir  and  (u*)s.  are  gov¬ 
erned  by  the  relative  6  dependences  of  these  contributions.  Since  these  dependences 
for  ((<«).?  and  lu„)r  are  the  same,  the  order-of-magnitude  comparison  of  the  two 
contributors  to  the  azimuthal  displacement  is  the  same  as  the  comparison  of  the 
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FlC.  2.  Dependence  of  computed  «*l«  murrains  on  source  depth,  which  u  given  *t  the  left  of  *j»-h 
theoretical  seismogram  See  caption  of  Figure  I  lor  further  details 


polar  and  azimuthal  -L  placement*  I  »»:n  loisionaJ  excitation,  i.e 
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The  results  of  our  direct  computation  of  (u#)r  and  lu,)r  are  shown  in  Figure  4. 
where  we  see  that  relative  to  the  azimuthal  displacements,  the  polar  seismogram  is 
a  straight  line.  In  fact,  the  numerical  results  show  that  the  ratio  of  maximum 
displacements  in  the  two  cases  is  about  HX*  to  1.  We  therefore  conclude  that  for 
normal  WWSSN  seismograms  which  contain  all  of  the  generated  energy  down  to  a 
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See  rapt  ion  uf  Figure  I  for  further  details. 


Fir..  4.  Companion  of  »r»muthnl  m  and  polar  »*•  r«*mi*»*nent«  ol  ifwpl^i  *m«ii  from  inweawl 
ei- nation  Tracer  l.tl  iUo*'ivte  actual  pouml  di«»4*«  e»«r*.t.  and  !ra«c*  m. >he  »•  t. 
inetrvment  described  in  the  v  apt  ion  of  Ficurr  I  and  tieuwoan  roll  •  filter  •  ■onp«ila'et«w»l  parar 
pine#  the  epicenter  i»  the  f«**K  of  Kamchatka  peiweul .  the  «utK>n  at  HNK.  *•»>.♦•  40*.  A  - 
km.  A  •  90*.  and  a  step-function  time  dependence  at  (he  focua. 


period  of  10  sec.  the  total  azimuthal  component  of  displacement  can  Se  computed 
b%  ■  onsidenng  only  the  torsional-wave  contribution 

Before  going  into  a  specific,  direct  comparison  of  theorencal  and  eKper1111ent.1l 
seismograms,  it  should  be  emphasized  that  we  have  not  made  any  s|wiial  attempt 
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to  fit  the  experimental  data  yet.  We  have  just  picked  the  hest  fit  to  the  experimental 
record  from  among  the  theoretical  traces  shown  in  Figures  I  to  3.  Even  so,  the 
feature-by-feature  agreement  is  striking  over  much  of  the  record.  This  is  shown  in 
Figure  5  The  first  point  to  note  is  the  success  of  the  theoretical  work  in  generating 
the  phases  ScS  and  sScS  that  are  reflected  from  the  core-mantle  boundary  This 
means  that  we  can  expect  this  kind  of  analysis  to  be  useful  over  the  entire  mantle 
down  to  the  core.  The  second  point  is  the  phase-hv-phnse  agreement  in  the  bo.lv 
waves  identified  in  the  figure,  and  the  decrease  in  the  quality  of  ihe  agi.  emo.t 
beyond  them.  According  to  the  recent  paper  by  Kauswl  ct  a l.  (1977).  the  properti-s 
ol  Sa  are  determined  by  the  strut  ’ ure  from  the  vicinity  of  the  "850-km '  discontinuity 
upward.  Therefore,  the  good  agreement  in  initial  body  waves,  and  the  relatively 
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•ttenuauon  model  (or  the  Pacific  that  is  *iv*n  by  Mitchel  i  th* 

poorer  agreement  in  surface  waves  implies  that  we  have  essentially  the  right 
structure  below  7t*t  km.  hut  that  we  need  some  ne-'d..  ■„  nUr  earth  model 

above  that  As  mentioned  inive.  ihe  structure  we  used  u  an  average  oceanic  model  ’ 
however,  the  Kamchatka- KNR  path  is  along  the  "old"  edge  of  the  Pacific  plate 
Recent  seismic  studies  (eg.  Leeds.  ,1  at  197s  Leeds  19791  mdicate  -hat  .he 
thickness  of  the  lid  grows  with  increasing  distance  from  the  rulges.  Tlierefore 
along  the  path,  the  lid  should  he  thicker  than  average  This  mav  be  one  source  of 
the  discrepancies  between  the  theoretical  and  experimental  iraces  in  Figure  5 
Another  mav  be  oceanic  crustal  thicknesses,  for  perhaps  luuu  km  north  ol  HMl 
that  may  exceed  30  km  (Furumoto  rr  at  ,  1973). 

To  give  some  indication  of  whether  modifications  in  the  structure  above  700  km 
can  be  expected  to  improve  the  comparison  in  Figure  S.  we  have  repeated  th- 
computations  with  three  modified  structure*.  The  structures  and  results  are  show  i 
in  Figure  8.  which  indicates  that  quite  a  bit  can  be  done  w.th  the  single  structure 
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before  making  the  obvious  generalization  of  the  technique  to  include 


interpretation - , 

lateral  heterogeneity.  Consideration  of  the  experimental  coda  is  beyond  the  scope 

of  this  letter^  -  — — - - — -  - - " 

Conclusions 

It  is  now  feasible  to  make  direct  comparisons  of  theoretical  and  experimental 
seisin  'grsms.  for  t.iisionnl-wave  records  obtainml  from  the  long-pei lod  insininu  ms 
of  the  WWSSN.  With  our  present  technique,  all  of  the  generated  en-rgy— at  penods 
above  10  sec— is  exhibited  on  the  calculated  lime  senes,  body  waves  as  well  as 
surface  waves.  Application  of  this  result  miRht  be  expected  to  provide  the  highest 
source  and  structural  resolution  that  has  yet  been  achieved. 

Our  latest  work  (Schwab  et  at..  1977)  with  spheroidal  waves  on  spherical,  gravi¬ 
tating  models  of  the  eanh.  shows  that  even  when  we  employ  highly  optimized 
versions  of  the  current  algorithms,  these  calculations  are  six  tunes  slower  than  the 
comparable  torsional-wave  computations.  We  conclude  from  'his  that  new  aieo- 
rithms  are  required  before  the  efficiency^  of  spheroidal-wave  computations  (includ¬ 
ing  gravity)  will  approach  a  level  justifying  computation  ot  the  associated  theoretical 
seismograms,  down  to  a  penod  as  low  as  10  sec. 
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